00001 #include "fixnodesel.h"
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011 double round_it(double number)
00012 {
00013 double fractpart,intpart;
00014 fractpart = modf (number,&intpart);
00015 if(fractpart >= 0.5){
00016 number=ceil(number);
00017 }
00018 else{
00019 number=floor(number);
00020 }
00021 return(number);
00022 }
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033
00034 bool isPrime(int number)
00035 {
00036 long i,max;
00037
00038 if (number < 2) return false;
00039 if (number == 2) return true;
00040 if (number % 2 == 0) return false;
00041
00042 max = (long) sqrt(number);
00043
00044 for (i = 3; i <= max; i += 2) {
00045 if (number % i == 0) {
00046 return false;
00047 }
00048 }
00049 return true;
00050 }
00051
00052
00053
00054
00055
00056
00057
00058
00059
00060
00061
00062
00063
00064 long compute_number_of_combination(long n, long k)
00065 {
00066 long i;
00067 long denom,nom,diff;
00068 long combination;
00069
00070 diff = n - k;
00071 nom = 1;
00072 for(i = diff+1; i < n+1; i++){
00073 nom *= i;
00074 }
00075 denom = 1;
00076 for(i = 2; i < k+1; i++){
00077 denom *= i;
00078 }
00079 combination = nom/denom;
00080
00081 return(combination);
00082 }
00083
00084
00085
00086
00087
00088
00089
00090
00091
00092
00093
00094
00095
00096
00097 fixnodesel::fixnodesel(int np,int mr,long nd,meshdescription meshd,long* domproces,char *nameproc, int nl,long mes)
00098 {
00099 mespr=mes;
00100
00101
00102 nproc=np;
00103
00104 myrank=mr;
00105
00106 ndom=nd;
00107
00108 domproc = domproces;
00109
00110
00111 strcpy (procName,nameproc);
00112 nameLength=nl;
00113
00114
00115 md = meshd;
00116
00117
00118 fictdom = NULL;
00119
00120 nodeidentif = NULL;
00121
00122 coarseidentif = NULL;
00123
00124 coarsenadjac = NULL;
00125
00126 coarseadjac = NULL;
00127
00128 curidentif=NULL;
00129
00130 surfidentif=NULL;
00131
00132 surfnod=NULL;
00133
00134 curadjac=NULL;
00135
00136 curnadjac=NULL;
00137
00138 curnod=NULL;
00139
00140 start=NULL;
00141
00142 end=NULL;
00143
00144 members=NULL;
00145
00146 nedges=NULL;
00147
00148 nmembers=NULL;
00149
00150 userdefnod=NULL;
00151
00152 nsurf = 0;
00153
00154 surfmembers = NULL;
00155
00156 nsurfmembers = NULL;
00157
00158 automember = NULL;
00159
00160 glinkdom = NULL;
00161
00162 glinknod = NULL;
00163
00164 loclinkdom = NULL;
00165
00166 coarseidentif = NULL;
00167
00168
00169 surfcenters = NULL;
00170
00171 nsurfcentres = NULL;
00172
00173 surfnodpoint = NULL;
00174
00175 surfnodmark = NULL;
00176
00177 nsurf = 0;
00178
00179 nuserdefnod = -1;
00180
00181 dim = -1;
00182
00183 minSelStrat = -1;
00184
00185 condfixing = nomethod;
00186
00187 typecondcur = notype;
00188
00189 nmembercur = 0;
00190
00191 typecondsurf = notype;
00192
00193 nmembersurf = 0;
00194
00195 methodcondcor = 0;
00196
00197 sizefictdom = -1.0;
00198
00199 nadjacboundnod = NULL;
00200
00201 adjacboundnod = NULL;
00202
00203
00204 tol = 1e-10;
00205
00206 ratio = NULL;
00207
00208 }
00209
00210
00211
00212
00213
00214
00215 fixnodesel::~fixnodesel()
00216 {
00217 long i,j;
00218
00219
00220 if(fictdom != NULL){
00221 delete []fictdom;
00222 }
00223 if( nodeidentif != NULL){
00224 delete []nodeidentif;
00225 }
00226 if( userdefnod != NULL){
00227 delete []userdefnod;
00228 }
00229
00230 if(myrank == 0){
00231 if(coarseadjac != NULL){
00232 for(i = 0; i < tnbn; i++){
00233 delete []coarseadjac[i];
00234 }
00235 }
00236 if(coarsenadjac != NULL){
00237 delete []coarsenadjac;
00238 }
00239
00240 if(coarseidentif != NULL){
00241 delete []coarseidentif;
00242 }
00243 if(surfidentif != NULL){
00244 delete []surfidentif;
00245 }
00246 if(curidentif != NULL){
00247 delete []curidentif;
00248 }
00249 if(curadjac != NULL){
00250 for(j = 0; j < ncurnod; j++){
00251 delete []curadjac[i];
00252 }
00253 delete []curadjac;
00254 }
00255 if(curnadjac != NULL){
00256 delete []curnadjac;
00257 }
00258 if(curnod != NULL){
00259 delete []curnod;
00260 }
00261 if( start != NULL){
00262 delete []start;
00263 }
00264 if( end != NULL){
00265 delete []end;
00266 }
00267 if( members != NULL){
00268 for(j = 0; j < ncurves; j++){
00269 delete []members[j];
00270 }
00271 delete []members;
00272 }
00273 if( nmembers != NULL){
00274 delete []nmembers;
00275 }
00276 if(nedges != NULL){
00277 delete []nedges;
00278 }
00279
00280 if(surfnod != NULL){
00281 delete []surfnod;
00282 }
00283 if(surfnodpoint != NULL){
00284
00285 }
00286 }
00287
00288 if(surfmembers != NULL){
00289 for(i = 0; i < nsurf; i++){
00290 delete surfmembers[i];
00291 }
00292 delete []surfmembers;
00293 }
00294 if(nsurfmembers != NULL){
00295 delete []nsurfmembers;
00296 }
00297
00298 if(glinkdom != NULL){
00299 for(i = 0; i < tnbn; i++){
00300 delete []glinkdom[i];
00301 }
00302 delete []glinkdom;
00303 }
00304 if(glinknod != NULL){
00305 for(i = 0; i < tnbn; i++){
00306 delete []glinknod[i];
00307 }
00308 delete []glinknod;
00309 }
00310
00311 if(loclinkdom != NULL){
00312 for(i = 0; i < nbn; i++){
00313 delete []loclinkdom[i];
00314 }
00315 delete []loclinkdom;
00316 }
00317
00318 if(automember != NULL){
00319 delete []automember;
00320 }
00321 if(surfcenters != NULL){
00322 for(i = 0; i < nsurf; i++){
00323 delete []surfcenters[i];
00324 }
00325 delete []surfcenters;
00326 }
00327 if(nsurfcentres != NULL){
00328 delete []nsurfcentres;
00329 }
00330
00331
00332 if(surfnodmark != NULL){
00333 for(i = 0; i < nsurf; i++){
00334 delete []surfnodmark[i];
00335 }
00336 delete []surfnodmark;
00337 }
00338
00339 delete []ratio;
00340 if(mespr == 1) fprintf(out,"\n\n\nFIXNODESEL CLASS - end\n");
00341 }
00342
00343
00344
00345
00346
00347
00348
00349
00350
00351
00352
00353 void fixnodesel::initiate(gtopology *topol,partop *part,FILE *outfile)
00354 {
00355
00356 out = outfile;
00357
00358 top = topol;
00359
00360 ptop = part;
00361
00362 nn = top->nn;
00363
00364 ne = top->ne;
00365
00366 if (myrank==0){
00367
00368 nbnd = ptop->nbnd;
00369
00370 multip = ptop->bmultip;
00371
00372
00373 cnbn = ptop->icnbnmas;
00374
00375 tnbn = ptop->tnbn;
00376 }
00377
00378 nodmultip = ptop->nodmultip;
00379
00380 nbn = ptop->nbn;
00381
00382 lnbn = ptop->lnbndom;
00383
00384 lgnbn = ptop->icnbn;
00385
00386 maxnbn = ptop->maxnbn;
00387
00388 nin = ptop->nin;
00389
00390 lnin= ptop->lnindom;
00391
00392 if(mespr == 1) fprintf(out,"\n\n\nFIXNODESEL CLASS - start\n");
00393
00394 }
00395
00396
00397
00398
00399
00400
00401
00402
00403
00404
00405
00406
00407 bool fixnodesel::check_ltg(long *ltg,long nnd,FILE *out)
00408 {
00409 long neg,pos,one,i;
00410 switch (md){
00411 case all_nodes:{
00412 return false;
00413 break;
00414 }
00415 case bound_nodes:{
00416 neg = 0;
00417 pos = 0;
00418 one = 0;
00419 for (i = 0; i < nnd; i++ ){
00420 if (ltg[i] <= -2){
00421 neg++;
00422 }
00423 if (ltg[i] > -1){
00424 pos++;
00425 }
00426 if (ltg[i] == -1){
00427 one++;
00428 }
00429 }
00430
00431
00432
00433
00434
00435 if ( neg > 0){
00436 if(mespr == 1) fprintf(out,"\n\n\nFixing nodes are found in ltg array\n");
00437 return false;
00438 }
00439
00440 else{
00441 if(mespr == 1) fprintf(out,"\n\n\nFixing nodes must be found by function fixing_detection\n");
00442 return true;
00443 }
00444 break;
00445 }
00446 case neg_bound_nodes:{
00447 return false;
00448 break;
00449 }
00450 default:{
00451
00452 }
00453 }
00454 fflush(out);
00455 }
00456
00457
00458
00459
00460
00461
00462
00463
00464 void fixnodesel::give_whole_dim()
00465 {
00466 long *elem;
00467 long i;
00468
00469 elem = new long[ne];
00470
00471 for(i = 0; i < ne; i++){
00472 elem[i] = top->give_whole_dim(i);
00473 }
00474 for(i = 1; i < ne; i++){
00475 if(elem[i] != elem[i-1]){
00476 par_print_err(myrank,"different spatial dimension is found in the mesh,function fixing_detection can not be used\n", __FILE__, __LINE__, __func__);
00477 }
00478 else{
00479 dim = elem[i];
00480 }
00481 }
00482
00483 delete []elem;
00484 }
00485
00486
00487
00488
00489
00490
00491
00492
00493
00494
00495
00496
00497
00498
00499
00500
00501
00502 double fixnodesel::compute_length_of_vector(long a,long b)
00503 {
00504 double u1,u2,u3,x1,x2,y1,y2,z1,z2,u;
00505
00506
00507 x1 = top->gnodes[a].x;
00508 y1 = top->gnodes[a].y;
00509 z1 = top->gnodes[a].z;
00510
00511 x2 = top->gnodes[b].x;
00512 y2 = top->gnodes[b].y;
00513 z2 = top->gnodes[b].z;
00514
00515 u1 = x2 - x1;
00516
00517 u2 = y2 - y1;
00518
00519 u3 = z2 - z1;
00520
00521 u = u1*u1+u2*u2+u3*u3;
00522 u=sqrt(u);
00523 return(u);
00524 }
00525
00526
00527
00528
00529
00530
00531
00532
00533
00534
00535
00536
00537
00538
00539
00540
00541 double fixnodesel::compute_angle_of_vector_a_and_b(long a,long b,long c)
00542 {
00543 double u1,u2,u3,v1,v2,v3,x1,x2,x3,y1,y2,y3,z1,z2,z3,scal,u,v,cosalpha,angle;
00544
00545
00546 x1 = top->gnodes[a].x;
00547 y1 = top->gnodes[a].y;
00548 z1 = top->gnodes[a].z;
00549
00550 x2 = top->gnodes[b].x;
00551 y2 = top->gnodes[b].y;
00552 z2 = top->gnodes[b].z;
00553
00554 u1 = x2 - x1;
00555
00556 u2 = y2 - y1;
00557
00558 u3 = z2 - z1;
00559
00560
00561 x3 = top->gnodes[c].x;
00562 y3 = top->gnodes[c].y;
00563 z3 = top->gnodes[c].z;
00564
00565 v1 = x3 - x1;
00566
00567 v2 = y3 - y1;
00568
00569 v3 = z3 - z1;
00570
00571
00572 scal = u1*v1+u2*v2+u3*v3;
00573
00574
00575 u = u1*u1+u2*u2+u3*u3;
00576
00577 u = sqrt(u);
00578
00579
00580 v = v1*v1+v2*v2+v3*v3;
00581
00582 v = sqrt(v);
00583
00584
00585 cosalpha = scal/(u*v);
00586
00587 if(cosalpha > 1.0){
00588 cosalpha = 1.0;
00589 }
00590 if(cosalpha < -1.0){
00591 cosalpha = -1.0;
00592 }
00593 angle = acos(cosalpha);
00594 return(angle);
00595
00596 }
00597
00598
00599
00600
00601
00602
00603
00604
00605
00606
00607
00608
00609
00610
00611
00612
00613 double fixnodesel::compute_area_of_triangle(long a,long b,long c)
00614 {
00615
00616 double x1,x2,x3,y1,y2,y3,z1,z2,z3,u1,u2,u3,v1,v2,v3,w1,w2,w3,w,area;
00617 double pr,dr;
00618
00619
00620
00621
00622
00623
00624
00625
00626
00627
00628
00629 x1 = top->gnodes[a].x;
00630 y1 = top->gnodes[a].y;
00631 z1 = top->gnodes[a].z;
00632
00633 x2 = top->gnodes[b].x;
00634 y2 = top->gnodes[b].y;
00635 z2 = top->gnodes[b].z;
00636
00637 u1 = x2 - x1;
00638
00639 u2 = y2 - y1;
00640
00641 u3 = z2 - z1;
00642
00643
00644 x3 = top->gnodes[c].x;
00645 y3 = top->gnodes[c].y;
00646 z3 = top->gnodes[c].z;
00647
00648 v1 = x3 - x1;
00649
00650 v2 = y3 - y1;
00651
00652 v3 = z3 - z1;
00653
00654
00655
00656 w1=u2*v3-u3*v2;
00657 w2=u3*v1-u1*v3;
00658 w3=u1*v2-u2*v1;
00659
00660 w = w1*w1+w2*w2+w3*w3;
00661
00662 area=w/2;
00663
00664
00665
00666
00667 return(area);
00668 }
00669
00670
00671
00672
00673
00674
00675 void fixnodesel::check_triangle()
00676 {
00677 long i,j,k;
00678 double *length,*angle;
00679 long v;
00680 long *vcomb,*ver,*pointver;
00681 long stop;
00682
00683
00684
00685 vcomb = new long[3];
00686
00687 v = 0 ;
00688 for(i = 0; i < nbn; i++){
00689 if(nodeidentif[i] == 3){
00690 v++;
00691 }
00692 }
00693 ver = new long[v];
00694 pointver = new long[v];
00695 j = 0 ;
00696 for(i = 0; i < nbn; i++){
00697 if(nodeidentif[i] == 3){
00698 ver[j] = lnbn[i];
00699 pointver[j] = i;
00700 j++;
00701 }
00702 }
00703
00704 length = new double[3];
00705 angle = new double[3];
00706 for(i = 0; i < v; i++){
00707 stop = 0;
00708 if(nodeidentif[pointver[i]] == 3){
00709 vcomb[0] = ver[i];
00710 for(j = i+1; j < v; j++){
00711 stop = 0;
00712 if(nodeidentif[pointver[j]] == 3){
00713 vcomb[1] = ver[j];
00714 length[2]= compute_length_of_vector(vcomb[0],vcomb[1]);
00715 for(k = j+1; k < v; k++){
00716 if(nodeidentif[pointver[k]] == 3){
00717 vcomb[2] = ver[k];
00718
00719 length[1] = compute_length_of_vector(vcomb[0],vcomb[2]);
00720 length[0] = compute_length_of_vector(vcomb[1],vcomb[2]);
00721 angle[0] = compute_angle_of_vector_a_and_b(vcomb[0],vcomb[1],vcomb[2]);
00722 length[0] = compute_length_of_vector(vcomb[1],vcomb[2]);
00723 angle[1] = compute_angle_of_vector_a_and_b(vcomb[1],vcomb[0],vcomb[2]);
00724 angle[2] = compute_angle_of_vector_a_and_b(vcomb[2],vcomb[0],vcomb[1]);
00725
00726
00727
00728
00729
00730
00731
00732
00733
00734
00735 if((angle[0] < 0.034906585 || angle[0] > 3.1066861) || (angle[1] < 0.034906585 || angle[1] > 3.1066861) || (angle[2] < 0.034906585 || angle[2] > 3.1066861)){
00736
00737 if(length[1] >= length[2]){
00738 nodeidentif[pointver[j]]=2;
00739 stop = 1;
00740
00741 }
00742 else{
00743 nodeidentif[pointver[k]]=2;
00744
00745 }
00746 }
00747
00748
00749
00750 if(angle[0] < (3.141592653589793238462643*5)/180 && angle[1] < (3.141592653589793238462643*5)/180){
00751 nodeidentif[pointver[k]]=2;
00752
00753 }
00754 if(angle[0] < (3.141592653589793238462643*5)/180 && angle[2] < (3.141592653589793238462643*5)/180 ){
00755 nodeidentif[pointver[j]]=2;
00756
00757
00758 }
00759 if(angle[1] < (3.141592653589793238462643*5)/180 && angle[2] < (3.141592653589793238462643*5)/180){
00760 nodeidentif[pointver[i]]=2;
00761
00762 }
00763
00764
00765 if(length[0]< 0.05*sizefictdom ){
00766 if(length[1] > 0.05*sizefictdom && length[1] > length[2]){
00767 nodeidentif[pointver[j]]=2;
00768
00769 stop = 1;
00770 }
00771 if(length[2] > 0.05*sizefictdom && length[2] > length[1]){
00772
00773 nodeidentif[pointver[k]]=2;
00774 }
00775 }
00776 if(length[1]< 0.05*sizefictdom ){
00777 if(length[0] > 0.05*sizefictdom && length[0] > length[2]){
00778
00779 nodeidentif[pointver[i]]=2;
00780 stop = 2;
00781 }
00782 if(length[2] > 0.05*sizefictdom && length[2] > length[0]){
00783 nodeidentif[pointver[k]]=2;
00784
00785 }
00786 }
00787 if(length[2]< 0.05*sizefictdom ){
00788 if(length[0] > 0.05*sizefictdom && length[0] > length[1]){
00789 nodeidentif[pointver[i]]=2;
00790
00791 stop = 2;
00792 }
00793 if(length[1] > 0.05*sizefictdom && length[1] > length[0]){
00794 nodeidentif[pointver[j]]=2;
00795
00796 stop = 1;
00797 }
00798 }
00799 if(stop == 1){
00800 break;
00801 }
00802 }
00803 }
00804 if(stop == 2){
00805 break;
00806 }
00807 }
00808 }
00809 }
00810 }
00811
00812 j = 0;
00813 if(mespr == 1){
00814 for(i = 0; i < nbn; i++){
00815
00816 if(nodeidentif[i] == 3){
00817 fprintf(out,"Node %ld is fixing nodes\n",lnbn[i]);
00818 j++;
00819 }
00820 }
00821 fprintf(out,"Number of fixing node is %ld, after check_triangle\n",j);
00822 }
00823 delete []vcomb;
00824 delete []ver;
00825 delete []pointver;
00826
00827 delete []length;
00828 delete []angle;
00829 }
00830
00831
00832
00833
00834
00835
00836
00837
00838
00839
00840
00841
00842
00843
00844
00845
00846 void fixnodesel::set_fictitious_subdomain()
00847 {
00848 long i;
00849 double miny,minx,minz,maxy,maxx,maxz,x,y,z;
00850
00851 if(fictdom != NULL){
00852 delete []fictdom;
00853 }
00854 fictdom = new double[6];
00855
00856
00857
00858
00859
00860
00861 minx = top->gnodes[0].x;
00862 miny = top->gnodes[0].y;
00863 minz = top->gnodes[0].z;
00864 maxx = top->gnodes[0].x;
00865 maxy = top->gnodes[0].y;
00866 maxz = top->gnodes[0].z;
00867 for(i = 0; i < nn; i++){
00868 x = top->gnodes[i].x;
00869 y = top->gnodes[i].y;
00870 z = top->gnodes[i].z;
00871
00872
00873 if(x <= minx){
00874 minx = x;
00875 fictdom[0] = x;
00876 }
00877
00878 if(y <= miny){
00879 miny = y;
00880 fictdom[1] = y;
00881 }
00882
00883 if(z <= minz){
00884 minz = z;
00885 fictdom[2] = z;
00886 }
00887
00888 if(x >= maxx){
00889 maxx = x;
00890 fictdom[3] = x;
00891 }
00892
00893 if(y >= maxy){
00894 maxy = y;
00895 fictdom[4] = y;
00896 }
00897
00898 if(z >= maxz){
00899 maxz = z;
00900 fictdom[5] = z;
00901 }
00902 }
00903 if(mespr == 1) fprintf(out,"fictdom: %lf %lf %lf %lf %lf %lf\n",fictdom[0],fictdom[1],fictdom[2],fictdom[3],fictdom[4],fictdom[5]);
00904 fflush(out);
00905 }
00906
00907
00908
00909
00910
00911
00912
00913
00914
00915
00916
00917
00918
00919
00920
00921
00922
00923
00924
00925
00926
00927
00928 void fixnodesel::compute_size_of_fictitious_subdomain()
00929 {
00930 double vetcx,vetcy,vetcz,sizevect;
00931
00932 vetcx = fictdom[3] - fictdom[0];
00933 vetcy = fictdom[4] - fictdom[1];
00934 vetcz = fictdom[5] - fictdom[2];
00935
00936 sizevect = vetcx*vetcx + vetcy*vetcy + vetcz*vetcz;
00937 sizefictdom = sqrt(sizevect);
00938 if(mespr == 1) fprintf(out,"size of fictdom %lf\n",sizefictdom);
00939 lx=fictdom[3]-fictdom[0];
00940 ly=fictdom[4]-fictdom[1];
00941 lz=fictdom[5]-fictdom[2];
00942
00943
00944 ratio = new double[3];
00945 if( ly != 0.0){
00946 ratio[0]=lx/ly;
00947 }
00948 else{
00949 ratio[0] = 0.0;
00950 }
00951 if( lz != 0.0){
00952 ratio[1]=lx/lz;
00953 }
00954 else{
00955 ratio[1] = 0.0;
00956 }
00957 if( lz != 0.0){
00958 ratio[2]=ly/ly;
00959 }
00960 else{
00961 ratio[2] = 0.0;
00962 }
00963 if(mespr == 1) fprintf(out,"size of fictdom: x side %lf, y side %lf,z side %lf\n",lx,ly,lz);
00964 if(mespr == 1) fprintf(out,"ratio: xy %lf, xz %lf,yz %lf\n",ratio[0],ratio[1],ratio[2]);
00965 fflush(out);
00966 }
00967
00968
00969
00970
00971
00972
00973
00974
00975
00976 void fixnodesel::compute_statistics_of_multiplicity()
00977 {
00978 long i,j,k,l;
00979 long buffsize;
00980 long *buff;
00981 MPI_Status stat;
00982
00983 if(mespr == 1) fprintf(out,"\n\n\nStatistics:\n");
00984
00985 maxmultip = 0;
00986 for(i = 0; i < nn; i++){
00987 if(nodmultip[i] > maxmultip){
00988 maxmultip = nodmultip[i];
00989 }
00990 }
00991
00992 if(myrank == 0){
00993
00994 nstatdom = new long [nproc];
00995 statdom = new long*[nproc];
00996
00997 j=domproc[0];
00998 nstatdom[j]=maxmultip;
00999 statdom[j] = new long[nstatdom[j]];
01000 for (i = 1; i < nproc; i++){
01001 MPI_Recv (&k,1,MPI_LONG,MPI_ANY_SOURCE,MPI_ANY_TAG,MPI_COMM_WORLD,&stat);
01002
01003
01004 j=domproc[stat.MPI_TAG];
01005 if (maxmultip < k){
01006 maxmultip = k;
01007 }
01008 nstatdom[j] = k;
01009 statdom[j] = new long[maxmultip];
01010 }
01011
01012 if(mespr == 1){
01013 fprintf(out,"\n\n\nMaximal multiplicity %ld\n",maxmultip);
01014 for (i = 0; i < nproc; i++){
01015 fprintf(out,"\n\nDomain %ld has maximal multiplicity %ld\n",i,nstatdom[i]);
01016 }
01017 }
01018
01019
01020 for (i=1;i<nproc;i++){
01021 MPI_Send (&maxmultip,1,MPI_LONG,i,myrank,MPI_COMM_WORLD);
01022 }
01023 }
01024 else{
01025 MPI_Send (&maxmultip,1,MPI_LONG,0,myrank,MPI_COMM_WORLD);
01026 MPI_Recv (&maxmultip,1,MPI_LONG,MPI_ANY_SOURCE,MPI_ANY_TAG,MPI_COMM_WORLD,&stat);
01027 }
01028 MPI_Barrier (MPI_COMM_WORLD);
01029
01030 buffsize = maxmultip;
01031 buff = new long[buffsize];
01032
01033 for(j = 0; j < maxmultip; j++){
01034 buff[j] = 0;
01035 for(i = 0; i < nn; i++){
01036 if(nodmultip[i] == j+1){
01037 buff[j]++;
01038 }
01039 }
01040
01041 }
01042
01043 if(myrank == 0){
01044 k=domproc[0];
01045 for(j = 0; j < maxmultip; j++){
01046 statdom[k][j] = buff[j];
01047 }
01048
01049 for (i = 1; i < nproc; i++){
01050 MPI_Recv (buff,buffsize,MPI_LONG,MPI_ANY_SOURCE,MPI_ANY_TAG,MPI_COMM_WORLD,&stat);
01051
01052
01053 k=domproc[stat.MPI_TAG];
01054 for(j = 0; j < maxmultip; j++){
01055 statdom[k][j] = buff[j];
01056 }
01057 }
01058
01059 if(mespr == 1){
01060 for (i = 0; i < nproc; i++){
01061 fprintf(out,"\n\nDomain %ld:\n",i);
01062 for(j= 0; j < maxmultip; j++){
01063 fprintf(out,"Number of nodes with multiplicity %ld is %ld\n",j+1,statdom[i][j]);
01064 }
01065 }
01066 }
01067
01068 long estMinCornod;
01069 l = 0;
01070 estMinCornod = nproc*3;
01071 if(mespr == 1) fprintf(out,"\n\nEstimate minimal number of fixing nodes in whole problem %ld\n",estMinCornod);
01072 for (i = 0; i < nproc; i++){
01073 k = 0;
01074 for(j= 1; j < maxmultip; j++){
01075 if(statdom[i][j] < estMinCornod){
01076 k++;
01077 }
01078 }
01079 if(k == 1){
01080 l++;
01081 }
01082 }
01083 if(l != 0){
01084 minSelStrat = 1;
01085 }
01086 else{
01087 minSelStrat = 2;
01088 }
01089 if(mespr == 1) fprintf(out,"\n\nselected strategy: %ld\n",minSelStrat);
01090 for (i = 1; i < nproc; i++){
01091 MPI_Send (&minSelStrat,1,MPI_LONG,i,myrank,MPI_COMM_WORLD);
01092 }
01093
01094 }
01095 else{
01096 MPI_Send (buff,buffsize,MPI_LONG,0,myrank,MPI_COMM_WORLD);
01097 MPI_Recv (&minSelStrat,1,MPI_LONG,MPI_ANY_SOURCE,MPI_ANY_TAG,MPI_COMM_WORLD,&stat);
01098 }
01099 MPI_Barrier (MPI_COMM_WORLD);
01100 delete []buff;
01101 }
01102
01103
01104
01105
01106
01107
01108
01109
01110
01111
01112
01113
01114
01115
01116 void fixnodesel::create_link_dom_nod()
01117 {
01118 long i,j,m,k,l,max;
01119 MPI_Status stat;
01120 long *buff,buffsize;
01121
01122
01123 if(myrank == 0){
01124
01125
01126 glinknod = new long*[tnbn];
01127 glinkdom = new long*[tnbn];
01128
01129 for(i = 0; i < tnbn; i++){
01130 glinknod[i] = new long[multip[i]];
01131 glinkdom[i] = new long[multip[i]];
01132 }
01133
01134 for(i = 0; i < tnbn; i++){
01135 m = 0;
01136 for(j = 0; j < nproc; j++){
01137 for(k = 0; k < nbnd[j]; k++){
01138 if(cnbn[j][k] == i){
01139 glinknod[i][m] = k;
01140 glinkdom[i][m] = j;
01141 m++;
01142 break;
01143 }
01144 }
01145 }
01146 }
01147 max = 0;
01148 for(i = 0; i < nproc; i++){
01149 m = 0;
01150 for(j = 0; j < tnbn; j++){
01151 for(k = 0; k < multip[j]; k++){
01152 if(glinkdom[j][k] == i){
01153 m += multip[j] - 1;
01154 }
01155 }
01156 }
01157 if(m > max){
01158 max = m;
01159 }
01160 }
01161
01162
01163 buffsize = max;
01164 for(i = 1; i < nproc; i++){
01165 MPI_Send (&buffsize,1,MPI_LONG,i,myrank,MPI_COMM_WORLD);
01166 }
01167
01168
01169
01170
01171
01172
01173
01174
01175
01176
01177
01178
01179
01180
01181
01182 }
01183
01184
01185
01186 else{
01187 MPI_Recv (&buffsize,1,MPI_LONG,MPI_ANY_SOURCE,MPI_ANY_TAG,MPI_COMM_WORLD,&stat);
01188 }
01189 MPI_Barrier (MPI_COMM_WORLD);
01190 max = buffsize;
01191
01192 buff = new long[buffsize];
01193
01194
01195 if(myrank == 0){
01196
01197 for(i = 1; i < nproc; i++){
01198 l = 0;
01199 for(j = 0; j < tnbn; j++){
01200 for(k =0; k < multip[j]; k++){
01201 if(glinkdom[j][k] == i){
01202 for(m = 0; m < k; m++){
01203 buff[l]=glinkdom[j][m];
01204 l++;
01205 }
01206 for(m = k+1; m < multip[j]; m++){
01207 buff[l]=glinkdom[j][m];
01208 l++;
01209 }
01210 break;
01211 }
01212 }
01213 }
01214 MPI_Send (buff,buffsize,MPI_LONG,i,myrank,MPI_COMM_WORLD);
01215 }
01216 l = 0;
01217 for(j = 0; j < tnbn; j++){
01218 for(k =0; k < multip[j]; k++){
01219 if(glinkdom[j][k] == 0){
01220 for(m = 0; m < k; m++){
01221 buff[l]=glinkdom[j][m];
01222 l++;
01223 }
01224 for(m = k+1; m < multip[j]; m++){
01225 buff[l]=glinkdom[j][m];
01226 l++;
01227 }
01228 break;
01229 }
01230 }
01231 }
01232 }
01233 else{
01234 MPI_Recv (buff,buffsize,MPI_LONG,MPI_ANY_SOURCE,MPI_ANY_TAG,MPI_COMM_WORLD,&stat);
01235 }
01236
01237 loclinkdom = new long*[nbn];
01238 l = 0;
01239 for(i = 0; i < nbn; i++){
01240 loclinkdom[i] = new long[nodmultip[lnbn[i]]-1];
01241 for(j = 0; j < nodmultip[lnbn[i]]-1; j++){
01242 loclinkdom[i][j]=buff[l];
01243 l++;
01244 }
01245 }
01246 delete []buff;
01247
01248
01249
01250
01251
01252
01253
01254
01255
01256 }
01257
01258
01259
01260
01261
01262
01263
01264
01265
01266
01267
01268
01269 void fixnodesel::create_subdom_graph_2d()
01270 {
01271 long i,j,k,l,m;
01272 long ned,nned,auxnedges;
01273 long *auxnned,*dupledges,*aux,*locbound;
01274 long **auxedgenodes;
01275
01276
01277 auxnedges = 0;
01278
01279 aux = new long[3];
01280 for(i = 0; i < ne; i++){
01281
01282 ned = top->give_ned(i);
01283
01284 nned = top->give_nned(i);
01285
01286 for(j = 0; j < ned; j++ ){
01287 top->give_edge_nodes(i,j,aux);
01288 m = 0;
01289 for(l = 0; l < nned; l++){
01290 if(nodmultip[aux[l]] > 1){
01291 m++;
01292 }
01293 }
01294 if(m > 1){
01295 auxnedges++;
01296 }
01297 }
01298 }
01299
01300
01301
01302
01303 auxnned = new long[auxnedges];
01304 auxedgenodes = new long*[auxnedges];
01305 dupledges = new long[auxnedges];
01306
01307
01308 k = 0;
01309 for(i = 0; i < ne; i++){
01310 ned = top->give_ned(i);
01311 nned = top->give_nned(i);
01312 for(j = 0; j < ned; j++ ){
01313 top->give_edge_nodes(i,j,aux);
01314 m = 0;
01315 for(l = 0; l < nned; l++){
01316 if(nodmultip[aux[l]] > 1){
01317 m++;
01318 }
01319 }
01320 if(m > 1){
01321 auxnned[k] = nned;
01322 auxedgenodes[k] = new long[nned];
01323 for(l = 0; l < nned; l++){
01324 auxedgenodes[k][l] = aux[l];
01325 }
01326
01327 if(nned == 2){
01328 if(auxedgenodes[k][0] > auxedgenodes[k][1]){
01329 l = auxedgenodes[k][1];
01330 auxedgenodes[k][1] = auxedgenodes[k][0];
01331 auxedgenodes[k][0] = l;
01332 }
01333 }
01334 if(nned == 3){
01335 if(auxedgenodes[k][0] > auxedgenodes[k][2]){
01336 l = auxedgenodes[k][2];
01337 auxedgenodes[k][2] = auxedgenodes[k][0];
01338 auxedgenodes[k][0] = l;
01339 }
01340 if(auxedgenodes[k][0] > auxedgenodes[k][1]){
01341 l = auxedgenodes[k][1];
01342 auxedgenodes[k][1] = auxedgenodes[k][0];
01343 auxedgenodes[k][0] = l;
01344 }
01345 if(auxedgenodes[k][1] > auxedgenodes[k][2]){
01346 l = auxedgenodes[k][2];
01347 auxedgenodes[k][2] = auxedgenodes[k][1];
01348 auxedgenodes[k][1] = l;
01349 }
01350 }
01351 dupledges[k] = 0;
01352 k++;
01353 }
01354 }
01355 }
01356 delete []aux;
01357
01358
01359
01360
01361
01362
01363
01364
01365
01366
01367
01368
01369
01370
01371 for(i = 0; i < auxnedges; i++){
01372 for(j = i+1; j < auxnedges; j++){
01373 if(auxedgenodes[i][0] == auxedgenodes[j][0]){
01374 if(auxnned[i] == 2 && auxnned[j] == 2){
01375 if(auxedgenodes[i][1] == auxedgenodes[j][1]){
01376 dupledges[j]++;
01377 dupledges[i]++;
01378 }
01379 }
01380 if(auxnned[i] == 3 && auxnned[i] == 3){
01381 if(auxedgenodes[i][1] == auxedgenodes[j][1]){
01382 if(auxedgenodes[i][2] == auxedgenodes[j][2]){
01383 dupledges[j]++;
01384 dupledges[i]++;
01385 }
01386 }
01387 }
01388 }
01389 }
01390 }
01391
01392
01393
01394
01395
01396
01397 locbound = new long[nn];
01398 m = 0;
01399 for(i = 0; i < nn; i++){
01400 if(nodmultip[i] > 1){
01401 locbound[i] = m;
01402 m++;
01403 }
01404 else{
01405 locbound[i] = -1;
01406 }
01407 }
01408
01409
01410
01411
01412
01413
01414
01415
01416
01417
01418
01419
01420
01421
01422
01423
01424
01425
01426
01427
01428
01429 nadjacboundnod = new long[nbn];
01430 for(i = 0; i < nbn; i++){
01431 nadjacboundnod[i] = 0;
01432 for(j = 0; j < auxnedges; j++){
01433 if(dupledges[j] == 0) {
01434 for(k = 0; k < auxnned[j]; k++){
01435 if(auxedgenodes[j][k] == lnbn[i]){
01436 auxedgenodes[j][k] = i;
01437 nadjacboundnod[i]++;
01438 }
01439 }
01440 }
01441 }
01442 }
01443
01444
01445
01446
01447
01448
01449
01450
01451
01452
01453
01454
01455
01456
01457 adjacboundnod = new long*[nbn];
01458 for(i = 0; i < nbn; i++){
01459 adjacboundnod[i] = new long[nadjacboundnod[i]];
01460 m = 0;
01461 for(j = 0; j < auxnedges; j++){
01462 if(dupledges[j] == 0) {
01463 if(auxnned[j] == 2){
01464 for(k = 0; k < auxnned[j]; k++){
01465 if(auxedgenodes[j][k] == i){
01466 if(k == 0){
01467 adjacboundnod[i][m] = auxedgenodes[j][1];
01468 m++;
01469 }
01470 if(k == 1){
01471 adjacboundnod[i][m] = auxedgenodes[j][0];
01472 m++;
01473 }
01474 }
01475 }
01476 }
01477 if(auxnned[j] == 3){
01478 for(k = 0; k < auxnned[j]; k++){
01479 if(auxedgenodes[j][k] == lnbn[i]){
01480 if(k == 0){
01481 adjacboundnod[i][m] = auxedgenodes[j][1];
01482 m++;
01483 }
01484 if(k == 1){
01485 adjacboundnod[i][m] = auxedgenodes[j][0];
01486 m++;
01487 adjacboundnod[i][m] = auxedgenodes[j][2];
01488 m++;
01489 }
01490 if(k == 2){
01491 adjacboundnod[i][m] = auxedgenodes[j][1];
01492 m++;
01493 }
01494 }
01495 }
01496 }
01497 }
01498 }
01499 }
01500
01501
01502
01503
01504
01505
01506
01507
01508
01509
01510
01511
01512 for(i = 0; i < nbn; i++){
01513 for(j = 0; j < nadjacboundnod[i]; j++){
01514 k = adjacboundnod[i][j];
01515 adjacboundnod[i][j] = lgnbn[k];
01516 }
01517 }
01518
01519
01520
01521 for(i = 0; i < nbn; i++){
01522 if(nadjacboundnod[i] == 2){
01523 if(adjacboundnod[i][0] > adjacboundnod[i][1]){
01524 j = adjacboundnod[i][1];
01525 adjacboundnod[i][1] = adjacboundnod[i][0];
01526 adjacboundnod[i][0] = j;
01527 }
01528 }
01529 }
01530
01531
01532
01533
01534
01535
01536
01537
01538
01539
01540
01541
01542 for(i = 0; i < auxnedges; i++){
01543 delete []auxedgenodes[i];
01544 }
01545 delete []auxnned;
01546 delete []auxedgenodes;
01547 delete []dupledges;
01548 delete []locbound;
01549
01550 }
01551
01552
01553
01554
01555
01556
01557
01558
01559
01560
01561 void fixnodesel::select_minimal_number_2d()
01562 {
01563 long i,auxnv;
01564
01565 if(myrank == 0){
01566
01567 if (coarseidentif != NULL)
01568 delete []coarseidentif;
01569 coarseidentif = new long [tnbn];
01570
01571 auxnv = 0;
01572
01573 for(i = 0; i < tnbn; i++){
01574
01575 if(multip[i] == 2 && coarsenadjac[i] == 1){
01576 coarseidentif[i] = 3;
01577 }
01578 if(multip[i] > 2 && coarsenadjac[i] > 2){
01579 coarseidentif[i] = 3;
01580 }
01581 if(multip[i] == 2 && coarsenadjac[i] == 2){
01582 coarseidentif[i] = 2;
01583 }
01584 if(coarseidentif[i] == 3){
01585 auxnv++;
01586 }
01587 }
01588 if(mespr == 1) fprintf(out,"Number of potencial fixing nodes in whole problem is %ld\n\n",auxnv);
01589 fflush(out);
01590
01591
01592 if(mespr == 1){
01593 for (i = 0; i < tnbn; i++){
01594 if(coarseidentif[i] == 3) fprintf(out,"\n node number %ld is fixing nodes node",i);
01595 }
01596 fflush(out);
01597 fprintf(out,"\n\n");
01598 }
01599 }
01600 }
01601
01602
01603
01604
01605
01606 void fixnodesel::check_minimal_number_2d()
01607 {
01608 long a,b,i,j,k,l,m,auxnv;
01609 long buffsize,*buff;
01610 MPI_Status stat;
01611 long *nvdom;
01612 long **stor;
01613
01614 buffsize = maxnbn+1;
01615 buff = new long[buffsize];
01616 if(myrank == 0){
01617 nvdom = new long[nproc];
01618 for(i = 0; i < nproc; i++){
01619 nvdom[i] = 0;
01620 }
01621
01622 for(i = 0; i < tnbn; i++){
01623 if(coarseidentif[i] == 3){
01624 for(j = 0; j < multip[i]; j++){
01625 k = glinkdom[i][j];
01626 nvdom[k]++;
01627 }
01628 }
01629 }
01630
01631
01632
01633
01634
01635 for(i = 1; i < nproc; i++){
01636 for(j = 0; j < nbnd[i]; j++){
01637 k = cnbn[i][j];
01638 buff[j] = coarseidentif[k];
01639 }
01640 buff[maxnbn] = nvdom[i];
01641 MPI_Send (buff,buffsize,MPI_LONG,i,myrank,MPI_COMM_WORLD);
01642 }
01643
01644 for(j = 0; j < nbnd[0]; j++){
01645 k = cnbn[0][j];
01646 buff[j] = coarseidentif[k];
01647 }
01648 buff[maxnbn] = nvdom[0];
01649 }
01650 else{
01651 MPI_Recv (buff,buffsize,MPI_LONG,MPI_ANY_SOURCE,MPI_ANY_TAG,MPI_COMM_WORLD,&stat);
01652 }
01653 MPI_Barrier (MPI_COMM_WORLD);
01654
01655 nodeidentif = new long[nbn];
01656
01657 auxnv = 0;
01658 for(i = 0; i < nbn; i++){
01659 nodeidentif[i] = buff[i];
01660 if(nodeidentif[i] == 3){
01661 auxnv++;
01662 if(mespr == 1) fprintf(out,"local number of fixng node: %ld\n\n",lnbn[i]);
01663 }
01664 }
01665
01666
01667 if(mespr == 1) fprintf(out,"Number of potential fixing nodes on subdomain is %ld\n\n",buff[maxnbn]);
01668 fflush(out);
01669
01670
01671 if(buff[maxnbn] < 2){
01672 auxnv = buff[maxnbn];
01673 select_fixing_nodes_geom_2d(auxnv);
01674 buff[maxnbn] = -2;
01675 }
01676 else{
01677 if(condfixing == nocondconer){
01678 check_triangle();
01679 auxnv = 0;
01680 for(i = 0; i < nbn; i++){
01681 if(nodeidentif[i] == 3){
01682 auxnv++;
01683 }
01684 }
01685 if(auxnv != buff[maxnbn]){
01686 buff[maxnbn] = -1*auxnv;
01687 }
01688 if(auxnv < 2){
01689 select_fixing_nodes_geom_2d(auxnv);
01690 buff[maxnbn] = -2;
01691 }
01692 }
01693 }
01694
01695 if(buff[maxnbn] < 0){
01696 for(i = 0; i < nbn; i++){
01697 buff[i] = nodeidentif[i];
01698 }
01699 }
01700 if(mespr == 1) fprintf(out,"Number of potential fixing nodes %d on subdomain %d\n\n",abs(buff[maxnbn]),myrank);
01701
01702
01703
01704 if(myrank == 0){
01705
01706 stor = new long*[nproc];
01707
01708 if (buff[maxnbn] < 0 ){
01709 nvdom[0] = -1*buff[maxnbn];
01710 stor[0] = new long [nbnd[0]];
01711 for(j = 0; j < nbnd[0]; j++){
01712 stor[0][j] = buff[j];
01713 }
01714 }
01715 else{
01716 stor[0]= new long [nbnd[0]];
01717 for(j = 0; j < nbnd[0]; j++){
01718 stor[0][j] = buff[j];
01719 }
01720 }
01721
01722
01723 for(i = 1; i < nproc; i++){
01724 MPI_Recv (buff,buffsize,MPI_LONG,MPI_ANY_SOURCE,MPI_ANY_TAG,MPI_COMM_WORLD,&stat);
01725 k=domproc[stat.MPI_TAG];
01726
01727 if (buff[maxnbn] < 0 ){
01728 nvdom[k] = -1*buff[maxnbn];
01729 stor[k] = new long [nbnd[k]];
01730 for(j = 0; j < nbnd[k]; j++){
01731 stor[k][j] = buff[j];
01732 }
01733 }
01734 else{
01735 stor[k] = new long [nbnd[k]];
01736 for(j = 0; j < nbnd[k]; j++){
01737 stor[k][j] = buff[j];
01738 }
01739 }
01740 }
01741
01742
01743 for(i = 0; i < tnbn; i++){
01744
01745 m = 0;
01746 l = 0;
01747 if (coarseidentif[i] == 3){
01748 a = 1;
01749 b = 0;
01750 }
01751 else{
01752 a = 0;
01753 b = 1;
01754 }
01755 for(j = 0; j < multip[i]; j++){
01756 if(coarseidentif[i] != stor[glinkdom[i][j]][glinknod[i][j]]){
01757 if(stor[glinkdom[i][j]][glinknod[i][j]] == 3){
01758 a++;
01759 }
01760 if(stor[glinkdom[i][j]][glinknod[i][j]] == 2){
01761 b++;
01762 }
01763 m++;
01764 }
01765 if(nvdom[glinkdom[i][j]] <= 2){
01766 l++;
01767 }
01768 }
01769 if( m != 0){
01770
01771 if(a >= b){
01772 if(mespr == 1) fprintf(out,"prepis %ld z 2 na 3\n",i);
01773 coarseidentif[i] = 3;
01774 for(j = 0; j < multip[i]; j++){
01775 if(stor[glinkdom[i][j]][glinknod[i][j]] == 2){
01776 nvdom[glinkdom[i][j]]++;
01777 }
01778 }
01779 }
01780 else{
01781 if(l == 0){
01782 if(mespr == 1) fprintf(out,"prepis %ld z 3 na 2\n",i);
01783 coarseidentif[i] = 2;
01784 for(j = 0; j < multip[i]; j++){
01785 if(stor[glinkdom[i][j]][glinknod[i][j]] == 3){
01786 nvdom[glinkdom[i][j]]--;
01787 }
01788 }
01789 }
01790 else{
01791
01792 }
01793 }
01794 }
01795 }
01796 if(mespr == 1){
01797 for(i = 0; i < nproc; i++){
01798 fprintf(out,"\n\n\n Number of fixing nodes on domain %ld is %ld\n",i,nvdom[i]);
01799 }
01800 }
01801 }
01802 else{
01803 MPI_Send (buff,buffsize,MPI_LONG,0,myrank,MPI_COMM_WORLD);
01804 }
01805 MPI_Barrier (MPI_COMM_WORLD);
01806
01807 delete []buff;
01808
01809
01810 if(myrank == 0){
01811 m = 0;
01812 for(i = 0; i < nproc; i++){
01813 nvdom[i] = 0;
01814 }
01815 for(i = 0; i < tnbn; i++){
01816 if(coarseidentif[i] == 3){
01817 m++;
01818 for(j = 0; j < multip[i]; j++){
01819 k = glinkdom[i][j];
01820 nvdom[k]++;
01821 }
01822 }
01823 }
01824 if(mespr == 1){
01825 fprintf(out,"\n\n\n Total number of fixing nodes in whole problem is %ld\n",m);
01826 for(i = 0; i < nproc; i++){
01827 fprintf(out,"\n\n\n Number of fixing nodes on domain %ld is %ld\n",i,nvdom[i]);
01828 }
01829 }
01830
01831 delete []nvdom;
01832 for(i = 0; i < nproc; i++){
01833 if(stor[i] != NULL){
01834 delete []stor[i];
01835 }
01836 }
01837 delete []stor;
01838 }
01839
01840 }
01841
01842
01843
01844
01845
01846 void fixnodesel::select_fixing_nodes_geom_2d(long auxnv)
01847 {
01848
01849 long i,j,max;
01850 long *ver;
01851 double maxradius,radius;
01852 long *vcomb;
01853
01854
01855 switch(auxnv){
01856 case 0:
01857 ver = new long[2];
01858
01859
01860
01861
01862
01863
01864
01865
01866
01867
01868
01869 srand((unsigned int) time(NULL));
01870 j = rand()%(nbn-1);
01871
01872 ver[0] = j;
01873
01874 vcomb = new long[2];
01875 vcomb[0] = lnbn[ver[0]];
01876 maxradius = 0.0;
01877 for(i = 0; i < nbn; i++){
01878 vcomb[1] = lnbn[i];
01879 radius = compute_length_of_vector(vcomb[0],vcomb[1]);
01880
01881 if(radius >= maxradius){
01882 maxradius = radius;
01883 j = i;
01884 }
01885 }
01886 ver[0] = j;
01887
01888
01889
01890 vcomb[0] = lnbn[ver[0]];
01891 maxradius = 0.0;
01892 for(i = 0; i < nbn; i++){
01893 vcomb[1] = lnbn[i];
01894 radius = compute_length_of_vector(vcomb[0],vcomb[1]);
01895 if(radius >= maxradius){
01896 j = i;
01897 maxradius = radius;
01898 }
01899 }
01900 ver[1] = j;
01901 nodeidentif[ver[0]]=3;
01902 nodeidentif[ver[1]]=3;
01903 vcomb[0] = lnbn[ver[0]];
01904 vcomb[1] = lnbn[ver[1]];
01905
01906
01907
01908 delete []ver;
01909 delete []vcomb;
01910
01911
01912
01913
01914
01915
01916
01917
01918 break;
01919 case 1:
01920 ver = new long[2];
01921 max = 0;
01922 for(i = 0; i < nbn; i++){
01923 if(nodeidentif[i] == 3){
01924 ver[0] = i;
01925 }
01926 }
01927
01928
01929
01930 vcomb = new long[2];
01931 vcomb[0] = lnbn[ver[0]];
01932 maxradius = 0.0;
01933 for(i = 0; i < nbn; i++){
01934 vcomb[1] = lnbn[i];
01935 radius = compute_length_of_vector(vcomb[0],vcomb[1]);
01936
01937 if(radius >= maxradius){
01938 maxradius = radius;
01939 j = i;
01940 }
01941 }
01942 ver[1] = j;
01943
01944
01945
01946
01947 nodeidentif[ver[1]]=3;
01948 vcomb[0] = lnbn[ver[0]];
01949 vcomb[1] = lnbn[ver[1]];
01950
01951 delete []ver;
01952 delete []vcomb;
01953
01954
01955
01956
01957
01958
01959 break;
01960 }
01961
01962 }
01963
01964
01965
01966
01967
01968
01969
01970
01971
01972
01973
01974
01975 void fixnodesel::create_subdom_graph_3d ()
01976 {
01977 long i,j,k,l,m;
01978 long ned,nned,auxnedges;
01979 long *auxnned,*dupledges,*aux;
01980 long **auxedgenodes;
01981 double ts,te;
01982
01983 ts = clock();
01984
01985 auxnedges = 0;
01986
01987 aux = new long[3];
01988 for(i = 0; i < ne; i++){
01989
01990 ned = top->give_ned(i);
01991 nned = top->give_nned(i);
01992 for(j = 0; j < ned; j++ ){
01993 top->give_edge_nodes(i,j,aux);
01994 m = 0;
01995 for(l = 0; l < nned; l++){
01996 if(nodmultip[aux[l]] > 1){
01997 m++;
01998 }
01999 }
02000 if(m > 1){
02001 auxnedges++;
02002 }
02003 }
02004 }
02005
02006 te = clock();
02007
02008
02009
02010 ts = clock();
02011 auxnned = new long[auxnedges];
02012 auxedgenodes = new long*[auxnedges];
02013 dupledges = new long[auxnedges];
02014 k = 0;
02015 for(i = 0; i < ne; i++){
02016 ned = top->give_ned(i);
02017 nned = top->give_nned(i);
02018 for(j = 0; j < ned; j++ ){
02019 top->give_edge_nodes(i,j,aux);
02020 m = 0;
02021 for(l = 0; l < nned; l++){
02022 if(nodmultip[aux[l]] > 1){
02023 m++;
02024 }
02025 }
02026 if(m > 1){
02027 auxnned[k] = nned;
02028 auxedgenodes[k] = new long[nned];
02029 for(l = 0; l < nned; l++){
02030 auxedgenodes[k][l] = aux[l];
02031 }
02032
02033 if(nned == 2){
02034 if(auxedgenodes[k][0] > auxedgenodes[k][1]){
02035 l = auxedgenodes[k][1];
02036 auxedgenodes[k][1] = auxedgenodes[k][0];
02037 auxedgenodes[k][0] = l;
02038 }
02039 }
02040 if(nned == 3){
02041 if(auxedgenodes[k][0] > auxedgenodes[k][2]){
02042 l = auxedgenodes[k][2];
02043 auxedgenodes[k][2] = auxedgenodes[k][0];
02044 auxedgenodes[k][0] = l;
02045 }
02046 if(auxedgenodes[k][0] > auxedgenodes[k][1]){
02047 l = auxedgenodes[k][1];
02048 auxedgenodes[k][1] = auxedgenodes[k][0];
02049 auxedgenodes[k][0] = l;
02050 }
02051 if(auxedgenodes[k][1] > auxedgenodes[k][2]){
02052 l = auxedgenodes[k][2];
02053 auxedgenodes[k][2] = auxedgenodes[k][1];
02054 auxedgenodes[k][1] = l;
02055 }
02056 }
02057 dupledges[k] = 0;
02058 k++;
02059 }
02060 }
02061 }
02062 delete []aux;
02063
02064
02065 te = clock();
02066
02067
02068
02069
02070
02071
02072
02073
02074
02075
02076
02077 ts = clock();
02078
02079 for(i = 0; i < auxnedges; i++){
02080 for(j = i+1; j < auxnedges; j++){
02081 if(auxedgenodes[i][0] == auxedgenodes[j][0]){
02082 if(auxnned[i] == 2 && auxnned[j] == 2){
02083 if(auxedgenodes[i][1] == auxedgenodes[j][1]){
02084 dupledges[j]++;
02085
02086 }
02087 }
02088 if(auxnned[i] == 3 && auxnned[i] == 3){
02089 if(auxedgenodes[i][1] == auxedgenodes[j][1]){
02090 if(auxedgenodes[i][2] == auxedgenodes[j][2]){
02091 dupledges[j]++;
02092
02093 }
02094 }
02095 }
02096 }
02097 }
02098 }
02099 te = clock();
02100
02101
02102
02103
02104
02105
02106 k = 0;
02107 for(i = 0; i < auxnedges; i++){
02108 if(dupledges[i] == 0){
02109 k++;
02110 }
02111 }
02112
02113 ts = clock();
02114
02115 nadjacboundnod = new long[nbn];
02116 for(i = 0; i < nbn; i++){
02117 nadjacboundnod[i] = 0;
02118 for(j = 0; j < auxnedges; j++){
02119 if(dupledges[j] == 0) {
02120 for(k = 0; k < auxnned[j]; k++){
02121 if(auxedgenodes[j][k] == lnbn[i]){
02122 auxedgenodes[j][k] = i;
02123 nadjacboundnod[i]++;
02124 }
02125 }
02126 }
02127 }
02128 }
02129 te = clock();
02130
02131
02132
02133
02134
02135 ts = clock();
02136
02137 adjacboundnod = new long*[nbn];
02138 for(i = 0; i < nbn; i++){
02139 adjacboundnod[i] = new long[nadjacboundnod[i]];
02140 m = 0;
02141 for(j = 0; j < auxnedges; j++){
02142 if(dupledges[j] == 0) {
02143 if(auxnned[j] == 2){
02144 for(k = 0; k < auxnned[j]; k++){
02145 if(auxedgenodes[j][k] == i){
02146 if(k == 0){
02147 adjacboundnod[i][m] = auxedgenodes[j][1];
02148 m++;
02149 }
02150 if(k == 1){
02151 adjacboundnod[i][m] = auxedgenodes[j][0];
02152 m++;
02153 }
02154 }
02155 }
02156 }
02157 if(auxnned[j] == 3){
02158 for(k = 0; k < auxnned[j]; k++){
02159 if(auxedgenodes[j][k] == lnbn[i]){
02160 if(k == 0){
02161 adjacboundnod[i][m] = auxedgenodes[j][1];
02162 m++;
02163 }
02164 if(k == 1){
02165 adjacboundnod[i][m] = auxedgenodes[j][0];
02166 m++;
02167 adjacboundnod[i][m] = auxedgenodes[j][2];
02168 m++;
02169 }
02170 if(k == 2){
02171 adjacboundnod[i][m] = auxedgenodes[j][1];
02172 m++;
02173 }
02174 }
02175 }
02176 }
02177 }
02178 }
02179 }
02180
02181
02182
02183
02184
02185
02186
02187
02188
02189
02190
02191 te = clock();
02192
02193
02194
02195 ts = clock();
02196
02197
02198 for(i = 0; i < nbn; i++){
02199 if(nadjacboundnod[i] == 2){
02200 if(adjacboundnod[i][0] > adjacboundnod[i][1]){
02201 j = adjacboundnod[i][1];
02202 adjacboundnod[i][1] = adjacboundnod[i][0];
02203 adjacboundnod[i][0] = j;
02204 }
02205 }
02206 else{
02207
02208 for(j = nadjacboundnod[i] - 1; j > 0; j--){
02209 for(k = 0; k < j; k++){
02210 if(adjacboundnod[i][k] > adjacboundnod[i][k+1]){
02211 m = adjacboundnod[i][k];
02212 adjacboundnod[i][k] = adjacboundnod[i][k+1];
02213 adjacboundnod[i][k+1] = m;
02214 }
02215 }
02216 }
02217 }
02218 }
02219 te = clock();
02220
02221
02222
02223
02224
02225
02226
02227
02228
02229
02230
02231 ts = clock();
02232
02233
02234 for(i = 0; i < nbn; i++){
02235 for(j = 0; j < nadjacboundnod[i]; j++){
02236 k = adjacboundnod[i][j];
02237 adjacboundnod[i][j] = lgnbn[k];
02238 }
02239 }
02240 te = clock();
02241
02242
02243 fflush(out);
02244
02245 for(i = 0; i < auxnedges; i++){
02246 delete []auxedgenodes[i];
02247 }
02248 delete []auxnned;
02249 delete []dupledges;
02250 delete []auxedgenodes;
02251 }
02252
02253
02254
02255
02256
02257
02258
02259
02260
02261
02262 void fixnodesel::select_minimal_number_3d()
02263 {
02264 long i,j,a,b,c,d;
02265 long auxnv;
02266
02267 if(myrank == 0){
02268
02269 if (coarseidentif != NULL)
02270 delete []coarseidentif;
02271
02272 coarseidentif = new long [tnbn];
02273 auxnv = 0;
02274 for(i = 0; i < tnbn; i++){
02275 coarseidentif[i] = 2;
02276 if(multip[i] > 2){
02277
02278 b = 0;
02279 c = 0;
02280 d = 0;
02281 for(j = 0; j < coarsenadjac[i]; j++){
02282 a = coarseadjac[i][j];
02283
02284 if(multip[a] < multip[i]){
02285 b++;
02286 }
02287 if(multip[a] == multip[i]){
02288 c++;
02289 }
02290 if(multip[a] > multip[i]){
02291 d++;
02292 }
02293 }
02294
02295
02296 if(c == 0 && d == 0 && b == coarsenadjac[i]){
02297 coarseidentif[i] = 3;
02298 }
02299
02300 if(c == 1 && d == 0 ){
02301 coarseidentif[i] = 3;
02302 }
02303 }
02304 if(coarseidentif[i] == 3){
02305 auxnv++;
02306 }
02307 }
02308 if(mespr == 1) fprintf(out,"Number of potencial fixing nodes in whole problem is %ld\n\n",auxnv);
02309 fflush(out);
02310 }
02311
02312 }
02313
02314
02315
02316
02317
02318
02319 void fixnodesel::select_fixing_nodes_on_master()
02320 {
02321 long i,j,k;
02322 long *auxnv,*buff;
02323 long acticmultip,stop,buffsize,estMinCornod;
02324 MPI_Status stat;
02325
02326
02327 buffsize = maxnbn;
02328 buff = new long[buffsize];
02329
02330 if(myrank == 0){
02331 estMinCornod = nproc*3;
02332
02333 auxnv = new long[nproc];
02334 for(i = 0; i < nproc; i++){
02335 auxnv[i]=0;
02336 }
02337
02338 coarseidentif = new long[tnbn];
02339
02340 for(i = 0; i < tnbn; i++){
02341 if(multip[i] == maxmultip){
02342 coarseidentif[i] = 3;
02343 for(j = 0; j < multip[i]; j++){
02344 auxnv[glinkdom[i][j]]++;
02345 }
02346 }
02347 else{
02348 coarseidentif[i] = 2;
02349 }
02350 }
02351
02352
02353
02354
02355
02356
02357 for(i = 0; i < nproc; i++){
02358 if(auxnv[i] < 3){
02359 stop = 0;
02360 acticmultip = maxmultip-1;
02361 while (stop != 1){
02362 if(statdom[i][acticmultip-1] != 0){
02363 for(j = 0; j < nbnd[i]; j++){
02364 if(multip[cnbn[i][j]] == acticmultip){
02365 coarseidentif[cnbn[i][j]] = 3;
02366 for(k = 0; k < multip[cnbn[i][j]]; k++){
02367 auxnv[glinkdom[cnbn[i][j]][k]]++;
02368 }
02369 }
02370 }
02371 if(auxnv[i] >= 3){
02372 stop = 1;
02373 }
02374 else{
02375 stop = 0;
02376 acticmultip--;
02377 }
02378 }
02379 else{
02380 acticmultip--;
02381 }
02382 }
02383 }
02384 }
02385
02386
02387
02388
02389
02390 j = 0;
02391 k = 0;
02392 for(i = 0; i < tnbn; i++){
02393 if(coarseidentif[i] == 3){
02394 j++;
02395 }
02396 else{
02397 k++;
02398 }
02399 }
02400 if(mespr == 1){
02401 fprintf(out,"\nNumber of fixing nodes in whole problem is %ld\n",j);
02402 fprintf(out,"\nNumber of boundary nodes in whole problem is %ld\n",k);
02403 }
02404 for(j = 1; j < nproc; j++){
02405 for(k = 0; k < nbnd[j]; k++){
02406 buff[k] = coarseidentif[cnbn[j][k]];
02407 }
02408 MPI_Send (buff,buffsize,MPI_LONG,j,myrank,MPI_COMM_WORLD);
02409 }
02410
02411 nodeidentif = new long[nbn];
02412 for(j = 0; j < nbnd[0]; j++){
02413 nodeidentif[j] = coarseidentif[cnbn[0][j]];;
02414 }
02415
02416 delete []auxnv;
02417 }
02418 else{
02419 MPI_Recv (buff,buffsize,MPI_LONG,MPI_ANY_SOURCE,MPI_ANY_TAG,MPI_COMM_WORLD,&stat);
02420
02421 nodeidentif = new long[nbn];
02422 for(i = 0; i < nbn; i++){
02423 nodeidentif[i] = buff[i];
02424 }
02425 }
02426 if(mespr == 1){
02427 for (i = 0; i < nbn; i++){
02428 if(nodeidentif[i] == 3) fprintf(out,"\n node number %ld is fixing node",lnbn[i]+1);
02429 }
02430 fflush(out);
02431 fprintf(out,"\n\n");
02432 }
02433 delete []buff;
02434 }
02435
02436
02437
02438
02439
02440 void fixnodesel::select_fixing_nodes_geom_3d(long auxnv)
02441 {
02442 long i,j,k,max;
02443 long *ver,first,second,third;
02444 double maxradius,radius,area,maxarea,cosalpha,diff;
02445 double *contrad;
02446 long *vcomb;
02447
02448
02449 if(auxnv == 0){
02450
02451 srand((unsigned int) time(NULL));
02452 j = rand()%(nbn-1);
02453
02454
02455
02456
02457
02458
02459
02460
02461
02462
02463
02464 vcomb = new long[3];
02465 vcomb[0] = lnbn[j];
02466 maxradius = 0.0;
02467 for(i = 0; i < nbn; i++){
02468 if(i != j){
02469 vcomb[1] = lnbn[i];
02470 radius = compute_length_of_vector(vcomb[0],vcomb[1]);
02471
02472 if(radius > maxradius){
02473 maxradius = radius;
02474 first = i;
02475 }
02476 }
02477 }
02478 nodeidentif[first]=3;
02479
02480
02481
02482 vcomb[0] = lnbn[first];
02483 maxradius = 0.0;
02484 for(i = 0; i < nbn; i++){
02485 if(i != first){
02486 vcomb[1] = lnbn[i];
02487 radius = compute_length_of_vector(vcomb[0],vcomb[1]);
02488
02489 if(radius > maxradius){
02490 second = i;
02491 maxradius = radius;
02492 }
02493 }
02494 }
02495 nodeidentif[second]=3;
02496 vcomb[0] = lnbn[first];
02497 vcomb[1] = lnbn[second];
02498
02499
02500
02501
02502 maxarea = 0.0;
02503 for(i = 0; i < nbn; i++){
02504 if(vcomb[0] != lnbn[i] && vcomb[1] != lnbn[i]){
02505 vcomb[2] = lnbn[i];
02506
02507 area = compute_area_of_triangle(vcomb[0],vcomb[1],vcomb[2]);
02508
02509 if(maxarea < area){
02510 maxarea = area;
02511 third = i;
02512 }
02513 }
02514 }
02515 nodeidentif[third] = 3;
02516
02517
02518 delete []vcomb;
02519
02520
02521
02522
02523
02524
02525
02526
02527 }
02528
02529 if(auxnv == 1){
02530 max = 0;
02531 for(i = 0; i < nbn; i++){
02532 if(nodeidentif[i] == 3){
02533 first = i;
02534 }
02535 }
02536
02537
02538
02539 vcomb = new long[3];
02540 vcomb[0] = lnbn[ver[0]];
02541 maxradius = 0.0;
02542 for(i = 0; i < nbn; i++){
02543 vcomb[1] = lnbn[i];
02544 radius = compute_length_of_vector(vcomb[0],vcomb[1]);
02545
02546 if(radius >= maxradius){
02547 maxradius = radius;
02548 second = i;
02549 }
02550 }
02551
02552
02553
02554
02555 nodeidentif[second]=3;
02556 vcomb[0] = lnbn[first];
02557 vcomb[1] = lnbn[second];
02558
02559 delete []ver;
02560
02561 maxarea = 0.0;
02562 for(i = 0; i < nbn; i++){
02563 vcomb[2] = lnbn[i];
02564 area = compute_area_of_triangle(vcomb[0],vcomb[1],vcomb[2]);
02565 if(maxarea < area){
02566 maxarea = area;
02567 }
02568 }
02569 for(i = 0; i < nbn; i++){
02570 vcomb[2] = lnbn[i];
02571 area = compute_area_of_triangle(vcomb[0],vcomb[1],vcomb[2]);
02572 if(maxarea < area){
02573 j = i;
02574 }
02575 }
02576 nodeidentif[j] = 3;
02577 delete []vcomb;
02578
02579
02580
02581
02582
02583
02584 }
02585 if(auxnv == 2){
02586 ver = new long[2];
02587 vcomb = new long[3];
02588 j = 0;
02589 k = 0;
02590 for(i = 0; i < nbn; i++){
02591 if(nodeidentif[i] == 3){
02592 ver[j] = lnbn[i];
02593 j++;
02594 }
02595 if(nodeidentif[i] == 2){
02596 k++;
02597 }
02598 }
02599 radius = compute_length_of_vector(ver[0],ver[1]);
02600
02601 maxradius = 0.0;
02602 j = 0;
02603 contrad = new double[k];
02604 vcomb[0]=ver[0];
02605 for(i = 0; i < nbn; i++){
02606 if(nodeidentif[i] == 2){
02607 vcomb[1] = lnbn[i];
02608 contrad[j] = compute_length_of_vector(vcomb[0],vcomb[1]);
02609
02610 if(maxradius < contrad[j]){
02611 maxradius = contrad[j];
02612 }
02613 j++;
02614 }
02615 }
02616
02617 j = 0;
02618 vcomb[0] = ver[1];
02619 for(i = 0; i < nbn; i++){
02620 if(nodeidentif[i] == 2){
02621
02622 if(maxradius < contrad[j]){
02623 vcomb[1] = lnbn[i];
02624 radius = compute_length_of_vector(vcomb[0],vcomb[1]);
02625
02626 if(radius > 0.05*sizefictdom){
02627 vcomb[2] = ver[0];
02628 cosalpha = compute_angle_of_vector_a_and_b(vcomb[0],vcomb[1],vcomb[2]);
02629 if(cosalpha > 0.034906585 || cosalpha < 3.1066861){
02630 nodeidentif[i] = 3;
02631
02632 }
02633 }
02634 }
02635 j++;
02636 }
02637 }
02638 delete []contrad;
02639 delete []ver;
02640 delete []vcomb;
02641 }
02642 auxnv = 3;
02643 }
02644
02645
02646
02647
02648
02649 void fixnodesel::check_minimal_number_3d()
02650 {
02651 long a,b,i,j,k,l,m,auxnv;
02652 long buffsize,*buff;
02653 MPI_Status stat;
02654 long *nvdom;
02655 long **stor;
02656
02657 buffsize = maxnbn+1;
02658 buff = new long[buffsize];
02659 if(myrank == 0){
02660 nvdom = new long[nproc];
02661 for(i = 0; i < nproc; i++){
02662 nvdom[i] = 0;
02663 }
02664
02665 for(i = 0; i < tnbn; i++){
02666 if(coarseidentif[i] == 3){
02667 for(j = 0; j < multip[i]; j++){
02668 k = glinkdom[i][j];
02669 nvdom[k]++;
02670 }
02671 }
02672 }
02673
02674
02675
02676
02677
02678 for(i = 1; i < nproc; i++){
02679 for(j = 0; j < nbnd[i]; j++){
02680 k = cnbn[i][j];
02681 buff[j] = coarseidentif[k];
02682 }
02683 buff[maxnbn] = nvdom[i];
02684 MPI_Send (buff,buffsize,MPI_LONG,i,myrank,MPI_COMM_WORLD);
02685 }
02686
02687 for(j = 0; j < nbnd[0]; j++){
02688 k = cnbn[0][j];
02689 buff[j] = coarseidentif[k];
02690 }
02691 buff[maxnbn] = nvdom[0];
02692 }
02693 else{
02694 MPI_Recv (buff,buffsize,MPI_LONG,MPI_ANY_SOURCE,MPI_ANY_TAG,MPI_COMM_WORLD,&stat);
02695 }
02696 MPI_Barrier (MPI_COMM_WORLD);
02697
02698 nodeidentif = new long[nbn];
02699
02700 for(i = 0; i < nbn; i++){
02701 nodeidentif[i] = buff[i];
02702 }
02703
02704 if(mespr == 1) fprintf(out,"Number of potential fixing nodes on subdomain is %ld\n\n",buff[maxnbn]);
02705 fflush(out);
02706
02707
02708 if(buff[maxnbn] < 3){
02709 auxnv = buff[maxnbn];
02710 select_fixing_nodes_geom_3d(auxnv);
02711 buff[maxnbn] = -3;
02712 }
02713 else{
02714 if(condfixing == nocondconer){
02715 check_triangle();
02716 auxnv = 0;
02717 for(i = 0; i < nbn; i++){
02718 if(nodeidentif[i] == 3){
02719 auxnv++;
02720 }
02721 }
02722 if(auxnv != buff[maxnbn]){
02723 buff[maxnbn] = -1*auxnv;
02724 }
02725 if(auxnv < 3){
02726 select_fixing_nodes_geom_3d(auxnv);
02727 buff[maxnbn] = -3;
02728 }
02729 }
02730 }
02731
02732 if(buff[maxnbn] < 0){
02733 for(i = 0; i < nbn; i++){
02734 buff[i] = nodeidentif[i];
02735 }
02736 }
02737
02738 if(mespr == 1) fprintf(out,"Number of potential fixing nodes %d on subdomain %d\n\n",abs(buff[maxnbn]),myrank);
02739
02740
02741
02742 if(myrank == 0){
02743
02744 stor = new long*[nproc];
02745
02746 if (buff[maxnbn] < 0 ){
02747 nvdom[0] = -1*buff[maxnbn];
02748 stor[0] = new long [nbnd[0]];
02749 for(j = 0; j < nbnd[0]; j++){
02750 stor[0][j] = buff[j];
02751 }
02752 }
02753 else{
02754 stor[0]= new long [nbnd[0]];
02755 for(j = 0; j < nbnd[0]; j++){
02756 stor[0][j] = buff[j];
02757 }
02758 }
02759
02760
02761 for(i = 1; i < nproc; i++){
02762 MPI_Recv (buff,buffsize,MPI_LONG,MPI_ANY_SOURCE,MPI_ANY_TAG,MPI_COMM_WORLD,&stat);
02763 k=domproc[stat.MPI_TAG];
02764
02765 if (buff[maxnbn] < 0 ){
02766 nvdom[k] = -1*buff[maxnbn];
02767 stor[k] = new long [nbnd[k]];
02768 for(j = 0; j < nbnd[k]; j++){
02769 stor[k][j] = buff[j];
02770 }
02771 }
02772 else{
02773 stor[k] = new long [nbnd[k]];
02774 for(j = 0; j < nbnd[k]; j++){
02775 stor[k][j] = buff[j];
02776 }
02777 }
02778 }
02779
02780
02781 for(i = 0; i < tnbn; i++){
02782
02783 m = 0;
02784 l = 0;
02785 if (coarseidentif[i] == 3){
02786 a = 1;
02787 b = 0;
02788 }
02789 else{
02790 a = 0;
02791 b = 1;
02792 }
02793 for(j = 0; j < multip[i]; j++){
02794 if(coarseidentif[i] != stor[glinkdom[i][j]][glinknod[i][j]]){
02795 if(stor[glinkdom[i][j]][glinknod[i][j]] == 3){
02796 a++;
02797 }
02798 if(stor[glinkdom[i][j]][glinknod[i][j]] == 2){
02799 b++;
02800 }
02801 m++;
02802 }
02803 if(nvdom[glinkdom[i][j]] <= 3){
02804 l++;
02805 }
02806 }
02807 if( m != 0){
02808
02809 if(a >= b){
02810
02811 coarseidentif[i] = 3;
02812 for(j = 0; j < multip[i]; j++){
02813 if(stor[glinkdom[i][j]][glinknod[i][j]] == 2){
02814 nvdom[glinkdom[i][j]]++;
02815 }
02816 }
02817 }
02818 else{
02819 if(l == 0){
02820
02821 coarseidentif[i] = 2;
02822 for(j = 0; j < multip[i]; j++){
02823 if(stor[glinkdom[i][j]][glinknod[i][j]] == 3){
02824 nvdom[glinkdom[i][j]]--;
02825 }
02826 }
02827 }
02828 else{
02829
02830 }
02831 }
02832 }
02833 }
02834 if(mespr == 1){
02835 for(i = 0; i < nproc; i++){
02836 fprintf(out,"\n\n\n Number of fixing nodes on domain %ld is %ld\n",i,nvdom[i]);
02837 }
02838 }
02839 }
02840 else{
02841 MPI_Send (buff,buffsize,MPI_LONG,0,myrank,MPI_COMM_WORLD);
02842 }
02843 MPI_Barrier (MPI_COMM_WORLD);
02844
02845 delete []buff;
02846
02847 if(mespr == 1){
02848 for(i = 0; i < nbn; i++){
02849 if(nodeidentif[i] == 3){
02850 fprintf(out,"local number of fixng node: %ld\n\n",lnbn[i]);
02851 }
02852 }
02853 }
02854
02855
02856
02857
02858 if(myrank == 0){
02859 m = 0;
02860 for(i = 0; i < nproc; i++){
02861 nvdom[i] = 0;
02862 }
02863 for(i = 0; i < tnbn; i++){
02864 if(coarseidentif[i] == 3){
02865 m++;
02866 for(j = 0; j < multip[i]; j++){
02867 k = glinkdom[i][j];
02868 nvdom[k]++;
02869 }
02870 }
02871 }
02872 if(mespr == 1){
02873 fprintf(out,"\n\n\n Total number of fixing nodes in whole problem is %ld\n",m);
02874 for(i = 0; i < nproc; i++){
02875 fprintf(out,"\n\n\n Number of fixing nodes on domain %ld is %ld\n",i,nvdom[i]);
02876 }
02877 }
02878
02879 delete []nvdom;
02880 for(i = 0; i < nproc; i++){
02881 if(stor[i] != NULL){
02882 delete []stor[i];
02883 }
02884 }
02885 delete []stor;
02886 }
02887 }
02888
02889
02890
02891
02892
02893
02894
02895
02896
02897
02898
02899
02900 void fixnodesel::create_master_graph()
02901 {
02902 long i,j,k,m,a,b;
02903 long buffsize,max,adr1,adr2;
02904 long *buff,**aux,*naux,**auxnadjac,*help;
02905 MPI_Status stat;
02906
02907 buffsize = nbn;
02908 for(i = 0; i < nbn; i++){
02909 buffsize += nadjacboundnod[i];
02910 }
02911
02912
02913
02914
02915 if (myrank==0){
02916
02917 max = buffsize;
02918
02919 for (i=1;i<nproc;i++){
02920 MPI_Recv (&k,1,MPI_LONG,MPI_ANY_SOURCE,MPI_ANY_TAG,MPI_COMM_WORLD,&stat);
02921
02922
02923 if (max<k){
02924 max=k;
02925 }
02926 }
02927
02928 buffsize = max;
02929
02930 for (i=1;i<nproc;i++){
02931 MPI_Send (&buffsize,1,MPI_LONG,i,myrank,MPI_COMM_WORLD);
02932 }
02933 }
02934
02935
02936
02937 else{
02938 MPI_Send (&buffsize,1,MPI_LONG,0,myrank,MPI_COMM_WORLD);
02939 MPI_Recv (&buffsize,1,MPI_LONG,MPI_ANY_SOURCE,MPI_ANY_TAG,MPI_COMM_WORLD,&stat);
02940 }
02941 MPI_Barrier (MPI_COMM_WORLD);
02942
02943
02944
02945 buff = new long[buffsize];
02946
02947 buff[0] = 0;
02948 for(i = 1; i < nbn; i++){
02949 buff[i] = buff[i-1]+nadjacboundnod[i-1];
02950 }
02951
02952 m = nbn;
02953 for(i = 0; i < nbn; i++){
02954 for(j = 0; j < nadjacboundnod[i]; j++){
02955 buff[m] = adjacboundnod[i][j];
02956 m++;
02957 }
02958 }
02959 for(i = m; i < buffsize; i++){
02960 buff[m] = -1;
02961 m++;
02962 }
02963
02964
02965
02966
02967
02968
02969
02970 for(i = 0; i < nbn; i++){
02971 delete []adjacboundnod[i];
02972 }
02973 delete []adjacboundnod;
02974 delete []nadjacboundnod;
02975
02976 if (myrank == 0){
02977 auxnadjac = new long*[nproc];
02978 aux = new long*[nproc];
02979 naux = new long[nproc];
02980
02981 k=domproc[0];
02982 auxnadjac[k] = new long[nbnd[k]+1];
02983 for(j = 0; j < nbnd[k]; j++){
02984 auxnadjac[k][j] = buff[j];
02985 }
02986 naux[k] = 0;
02987 m = 0;
02988 for(j = nbnd[k]; j < buffsize; j++){
02989 if(buff[j] != -1){
02990 naux[k]++;
02991 }
02992 else{
02993 m++;
02994 }
02995 if(m == 2){
02996 break;
02997 }
02998 }
02999
03000 aux[k] = new long[naux[k]];
03001 auxnadjac[k][nbnd[k]] = naux[k];
03002 for(j = 0; j < naux[k]; j++){
03003 aux[k][j] = buff[j+nbnd[k]];
03004 }
03005
03006
03007
03008 for (i=1;i<nproc;i++){
03009 MPI_Recv (buff,buffsize,MPI_LONG,MPI_ANY_SOURCE,MPI_ANY_TAG,MPI_COMM_WORLD,&stat);
03010 k=domproc[stat.MPI_TAG];
03011 auxnadjac[k] = new long[nbnd[k]+1];
03012 for(j = 0; j < nbnd[k]; j++){
03013 auxnadjac[k][j] = buff[j];
03014 }
03015 naux[k] = 0;
03016 m = 0;
03017 for(j = nbnd[k]; j < buffsize; j++){
03018 if(buff[j] != -1){
03019 naux[k]++;
03020 }
03021 else{
03022 m++;
03023 }
03024 if(m == 2){
03025 break;
03026 }
03027 }
03028
03029 aux[k] = new long[naux[k]];
03030 auxnadjac[k][nbnd[k]] = naux[k];
03031 for(j = 0; j < naux[k]; j++){
03032 aux[k][j] = buff[j+nbnd[k]];
03033 }
03034 }
03035
03036
03037
03038
03039
03040
03041
03042
03043
03044
03045
03046
03047
03048
03049
03050
03051
03052
03053
03054
03055
03056
03057 max = 0;
03058 for(i = 0; i < tnbn; i++){
03059 b = 0;
03060 for(j = 0; j < multip[i]; j++){
03061 adr1=auxnadjac[glinkdom[i][j]][glinknod[i][j]];
03062 adr2=auxnadjac[glinkdom[i][j]][glinknod[i][j]+1];
03063 a = adr2-adr1;
03064 b+=a;
03065 }
03066 if(max < b){
03067 max = b;
03068 }
03069 }
03070
03071
03072 coarsenadjac = new long[tnbn];
03073 coarseadjac = new long*[tnbn];
03074 help = new long[max];
03075 for(i = 0; i < tnbn; i++){
03076
03077 b = 0;
03078 for(j = 0; j < multip[i]; j++){
03079 adr1=auxnadjac[glinkdom[i][j]][glinknod[i][j]];
03080 adr2=auxnadjac[glinkdom[i][j]][glinknod[i][j]+1];
03081 for(k = adr1; k < adr2; k++){
03082 help[b] = aux[glinkdom[i][j]][k];
03083 b++;
03084 }
03085 }
03086
03087
03088
03089
03090
03091
03092 if(b == 2){
03093 if(help[0] != help[1]){
03094 par_print_err(myrank,"problem with node %ld\n\n", __FILE__, __LINE__, __func__);
03095 }
03096 else{
03097 coarsenadjac[i] = 1;
03098 coarseadjac[i] = new long[coarsenadjac[i]];
03099 coarseadjac[i][0] = help[0];
03100
03101 }
03102 }
03103 else{
03104
03105 for(j = b - 1; j > 0; j--){
03106 for(k = 0; k < j; k++){
03107 if(help[k] > help[k+1]){
03108 m = help[k];
03109 help[k] = help[k+1];
03110 help[k+1] = m;
03111 }
03112 }
03113 }
03114
03115
03116
03117
03118
03119
03120 coarsenadjac[i] = 1;
03121 for(j = 1; j < b; j++){
03122 if(help[j] != help[j-1]){
03123 coarsenadjac[i]++;
03124 }
03125 }
03126
03127 coarseadjac[i] = new long[coarsenadjac[i]];
03128 coarseadjac[i][0] = help[0];
03129 m = 1;
03130 for(j = 1; j < b; j++){
03131 if(help[j] != help[j-1]){
03132 coarseadjac[i][m] = help[j];
03133 m++;
03134 }
03135 }
03136 }
03137 }
03138 delete []help;
03139
03140
03141
03142
03143
03144
03145
03146
03147
03148 for(i = 0; i < nproc; i++){
03149 delete []aux[i];
03150 delete []auxnadjac[i];
03151 }
03152 delete []aux;
03153 delete []naux;
03154 delete []auxnadjac;
03155 }
03156 else{
03157 MPI_Send (buff,buffsize,MPI_LONG,0,myrank,MPI_COMM_WORLD);
03158 }
03159 MPI_Barrier (MPI_COMM_WORLD);
03160 delete [] buff;
03161 }
03162
03163
03164
03165
03166
03167
03168
03169
03170
03171
03172
03173 void fixnodesel::assemble_ltg_with_fixings(long *ltg)
03174 {
03175 long i,k;
03176 long buffsize,c,b;
03177 long *buff,*coarsenumb;
03178 MPI_Status stat;
03179
03180
03181
03182 buffsize = maxnbn;
03183 buff = new long[buffsize];
03184
03185
03186 if(myrank==0){
03187 coarsenumb = new long[tnbn];
03188
03189 c = 1;
03190 b = 1;
03191 for(i = 0; i < tnbn; i++){
03192
03193 if(coarseidentif[i] == 3){
03194 coarsenumb[i] = c*(-1);
03195 c++;
03196
03197 }
03198 if(coarseidentif[i] == 2){
03199 coarsenumb[i] = b;
03200 b++;
03201
03202 }
03203 }
03204 if(mespr == 1){
03205 fprintf(out,"\n\n\n Number of fixing nodes in whole problem is %ld\n",c-1);
03206 fprintf(out,"\n\n\n Number of boundary nodes in whole problem is %ld\n",b-1);
03207 }
03208 if((c-1)+(b-1) != tnbn){
03209 par_print_err(myrank,"Problem: the number of fixing and remainning nodes is not equal to total number of boundary nodes\n", __FILE__, __LINE__, __func__);
03210 if(mespr == 1) fprintf (out,"Problem: the number of fixing and remainning nodes is not equal to total number of boundary nodes (file %s, line %d).\n",__FILE__,__LINE__);
03211
03212 }
03213
03214 for (i = 1;i < nproc; i++){
03215 for(k = 0; k < nbnd[i]; k++){
03216 buff[k] = coarsenumb[cnbn[i][k]];
03217 }
03218 MPI_Send (buff,buffsize,MPI_LONG,i,myrank,MPI_COMM_WORLD);
03219 }
03220
03221
03222 for(k = 0; k < nbnd[0]; k++){
03223 buff[k] = coarsenumb[cnbn[0][k]];
03224 }
03225 delete []coarsenumb;
03226 }
03227
03228
03229 else{
03230 MPI_Recv (buff,buffsize,MPI_LONG,MPI_ANY_SOURCE,MPI_ANY_TAG,MPI_COMM_WORLD,&stat);
03231 }
03232 MPI_Barrier (MPI_COMM_WORLD);
03233
03234
03235
03236
03237
03238
03239
03240 for(i = 0; i < nbn; i++){
03241 ltg[lnbn[i]]=buff[i]-1;
03242 }
03243 nin = nn - nbn;
03244 for(i = 0; i < nin; i++){
03245 ltg[lnin[i]]=-1;
03246 }
03247
03248
03249
03250
03251
03252 fflush(out);
03253
03254 delete []buff;
03255 }
03256
03257
03258
03259
03260
03261
03262
03263
03264
03265
03266
03267
03268
03269 void fixnodesel::select_centers_of_surfaces()
03270 {
03271 long i,j,k,l,max;
03272 long buffsize,*buff,*aux;
03273 MPI_Status stat;
03274 double t1,t2,tc;
03275
03276 long n;
03277 long **centres;
03278 long *ncentres;
03279 double x,y,z,cx,cy,cz,nc;
03280 double min;
03281 double *lengths;
03282 long **auxnodes;
03283 long subdomnsurf;
03284 long *subdomnsurfmem;
03285 long **subdomsurfmem;
03286
03287
03288 long maxncentres;
03289 maxncentres = 5;
03290
03291 buffsize = maxnbn;
03292 buff = new long[buffsize];
03293
03294
03295 if(myrank == 0){
03296 for(i = 1; i < nproc; i++){
03297 for(j = 0; j < nbnd[i]; j++){
03298 buff[j] = surfnod[cnbn[i][j]];
03299 }
03300 MPI_Send (buff,buffsize,MPI_LONG,i,myrank,MPI_COMM_WORLD);
03301
03302 }
03303 for(j = 0; j < nbnd[0]; j++){
03304 buff[j] = surfnod[cnbn[0][j]];
03305 }
03306 }
03307 else{
03308 MPI_Recv (buff,buffsize,MPI_LONG,MPI_ANY_SOURCE,MPI_ANY_TAG,MPI_COMM_WORLD,&stat);
03309 }
03310 MPI_Barrier (MPI_COMM_WORLD);
03311
03312 max = 0;
03313 for(i = 0; i < nbn; i++){
03314 if(max < buff[i]){
03315 max = buff[i];
03316 }
03317 }
03318
03319
03320 max++;
03321 aux = new long[max];
03322 for(i = 0; i < max; i++){
03323 aux[i] = 0;
03324 }
03325
03326 for(i = 0; i < nbn; i++){
03327 if(buff[i] != -1){
03328 j = buff[i];
03329 aux[j]++;
03330 }
03331 }
03332
03333
03334 subdomnsurf = 0;
03335 for(i = 0; i < max; i++){
03336 if(aux[i] != 0){
03337 subdomnsurf++;
03338 }
03339 }
03340
03341 subdomnsurfmem = new long[subdomnsurf];
03342 subdomsurfmem = new long*[subdomnsurf];
03343 for(i = 0; i < subdomnsurf; i++){
03344 subdomnsurfmem[i] = 0;
03345 }
03346
03347 j = 0;
03348 for(i = 0; i < max; i++){
03349 if(aux[i] != 0){
03350 subdomsurfmem[j] = new long[aux[i]];
03351 aux[i] = j;
03352 j++;
03353 }
03354 }
03355
03356 for(i = 0; i < nbn; i++){
03357 if(buff[i] != -1){
03358 j = buff[i];
03359 k = aux[j];
03360 subdomsurfmem[k][subdomnsurfmem[k]]=i;
03361 subdomnsurfmem[k]++;
03362 buff[i] = -2;
03363 }
03364 }
03365 delete[]aux;
03366
03367
03368
03369
03370
03371
03372
03373
03374
03375
03376
03377 max = 0;
03378 for(i = 0; i < subdomnsurf; i++){
03379 if(max < subdomnsurfmem[i]){
03380 max = subdomnsurfmem[i];
03381 }
03382 }
03383
03384
03385 lengths = new double[max];
03386 ncentres = new long[subdomnsurf];
03387 centres = new long*[subdomnsurf];
03388 realcg = new double*[subdomnsurf];
03389 for(i = 0; i < subdomnsurf; i++){
03390 realcg[i] = new double[3];
03391 tc = 0.0;
03392 t1= clock();
03393 x = 0.0;
03394 y = 0.0;
03395 z = 0.0;
03396 nc = 0.0;
03397 for(j = 0; j < subdomnsurfmem[i]; j++){
03398 x += top->gnodes[lnbn[subdomsurfmem[i][j]]].x;
03399 y += top->gnodes[lnbn[subdomsurfmem[i][j]]].y;
03400 z += top->gnodes[lnbn[subdomsurfmem[i][j]]].z;
03401 nc += 1.00;
03402 }
03403
03404
03405
03406
03407 cx=x/nc;
03408 cy=y/nc;
03409 cz=z/nc;
03410 realcg[i][0]=cx;
03411 realcg[i][1]=cy;
03412 realcg[i][2]=cz;
03413
03414
03415
03416
03417 min=DBL_MAX;
03418
03419
03420
03421
03422
03423 for(j = 0; j < subdomnsurfmem[i]; j++){
03424 x = top->gnodes[lnbn[subdomsurfmem[i][j]]].x;
03425 y = top->gnodes[lnbn[subdomsurfmem[i][j]]].y;
03426 z = top->gnodes[lnbn[subdomsurfmem[i][j]]].z;
03427 lengths[j]=(cx-x)*(cx-x)+(cy-y)*(cy-y)+(cz-z)*(cz-z);
03428
03429
03430 if(lengths[j] < min){
03431 min = lengths[j];
03432
03433 }
03434 }
03435
03436
03437
03438 double diff;
03439 ncentres[i] = 0;
03440 for(j = 0; j < subdomnsurfmem[i]; j++){
03441 diff=fabs(min-lengths[j]);
03442 if(diff < tol){
03443 ncentres[i]++;
03444 }
03445 }
03446
03447
03448 centres[i] = new long[ncentres[i]];
03449 n = 0;
03450
03451 for(j = 0; j < subdomnsurfmem[i]; j++){
03452 diff=fabs(min-lengths[j]);
03453 if(diff < tol){
03454
03455 centres[i][n] = subdomsurfmem[i][j];
03456 n++;
03457 }
03458 }
03459 t2= clock();
03460
03461 tc += ((t2-t1)/(double)CLOCKS_PER_SEC);
03462
03463
03464
03465
03466 }
03467 delete []lengths;
03468
03469
03470
03471
03472
03473
03474
03475 if(mespr == 1){
03476 for(i = 0; i < subdomnsurf; i++){
03477 fprintf(out,"Surface %ld has %ld centres:\n",i,ncentres[i]);
03478 for(j = 0; j < ncentres[i]; j++){
03479 fprintf(out," %ld",centres[i][j]);
03480 }
03481 fprintf(out,"\n");
03482 }
03483 }
03484
03485
03486 for(i = 0; i < subdomnsurf; i++){
03487 for(j = 0; j < ncentres[i]; j++){
03488 k = centres[i][j];
03489 buff[k] = 1;
03490 }
03491 }
03492
03493
03494
03495
03496
03497 if(myrank == 0){
03498 auxnodes = new long*[nproc];
03499 auxnodes[0] = new long[nbnd[0]];
03500 for(j = 0; j < nbnd[0]; j++){
03501 auxnodes[0][j] = buff[j];
03502 }
03503 for(i = 1; i < nproc; i++){
03504 MPI_Recv (buff,buffsize,MPI_LONG,MPI_ANY_SOURCE,MPI_ANY_TAG,MPI_COMM_WORLD,&stat);
03505 k = domproc[stat.MPI_TAG];
03506 auxnodes[k] = new long[nbnd[k]];
03507 for(j = 0; j < nbnd[k]; j++){
03508 auxnodes[k][j] = buff[j];
03509 }
03510 }
03511 }
03512 else{
03513 MPI_Send (buff,buffsize,MPI_LONG,0,myrank,MPI_COMM_WORLD);
03514 }
03515 MPI_Barrier (MPI_COMM_WORLD);
03516 delete []buff;
03517
03518 if(myrank == 0){
03519
03520
03521
03522
03523
03524
03525
03526 nsurfcentres = new long[nsurf];
03527 surfcenters = new long*[nsurf];
03528 for(i = 0; i < nsurf; i++){
03529 nsurfcentres[i] = 0;
03530 }
03531
03532 for(i = 0; i < nproc; i++){
03533 for(j = 0; j < nbnd[i]; j++){
03534 if(auxnodes[i][j] == 1){
03535 k = cnbn[i][j];
03536 l = surfnod[k];
03537 if(glinkdom[k][1] <= i){
03538 nsurfcentres[l]++;
03539 }
03540 }
03541 }
03542 }
03543
03544
03545 for(i = 0; i < nsurf; i++){
03546
03547 surfcenters[i] = new long[nsurfcentres[i]];
03548 nsurfcentres[i] = 0;
03549 }
03550
03551 for(i = 0; i < nproc; i++){
03552 for(j = 0; j < nbnd[i]; j++){
03553 if(auxnodes[i][j] == 1){
03554 k = cnbn[i][j];
03555 l = surfnod[k];
03556 if(glinkdom[k][1] <= i){
03557 surfcenters[l][nsurfcentres[l]]=k;
03558 nsurfcentres[l]++;
03559 }
03560 }
03561 }
03562 }
03563 if(mespr == 1){
03564 for(i = 0; i < nsurf; i++){
03565 fprintf(out,"Surface %ld has %ld centres:\n",i,nsurfcentres[i]);
03566 for(j = 0; j < nsurfcentres[i]; j++){
03567 fprintf(out," %ld",surfcenters[i][j]);
03568 }
03569 fprintf(out,"\n");
03570 }
03571 }
03572 for(i = 0; i < nproc; i++){
03573 delete []auxnodes[i];
03574 }
03575 delete []auxnodes;
03576 }
03577
03578 for(i = 0; i < subdomnsurf; i++){
03579 delete []centres[i];
03580 }
03581 delete []centres;
03582 delete []ncentres;
03583
03584 }
03585
03586
03587
03588
03589
03590
03591
03592
03593
03594
03595
03596 void fixnodesel::mark_surfnodes()
03597 {
03598 long i,j,k,l,m,n,a,b;
03599
03600 surfnodmark = new long*[nsurf];
03601 for(i = 0; i < nsurf; i++){
03602
03603 surfnodmark[i] = new long[nsurfmembers[i]];
03604
03605 for(j = 0; j < nsurfmembers[i]; j++){
03606 surfnodmark[i][j] = -1;
03607 }
03608
03609 for(j = 0; j < nsurfcentres[i]; j++){
03610 k = surfcenters[i][j];
03611 for(l = 0; l < nsurfmembers[i]; l++){
03612 if(surfmembers[i][l] == k){
03613 surfnodmark[i][l] = 0;
03614 break;
03615 }
03616 }
03617 }
03618
03619 m = 0;
03620 n = 0;
03621 for(m = 0; m < nsurfmembers[i]; m++){
03622 for(j = 0; j < nsurfmembers[i]; j++){
03623 if(surfnodmark[i][j] == m){
03624 a = surfmembers[i][j];
03625 b = coarsenadjac[a];
03626
03627 for(k = 0; k < b; k++){
03628 l = coarseadjac[a][k];
03629
03630 if(surfnod[l] == i && surfnodmark[i][surfnodpoint[l]] == -1){
03631 surfnodmark[i][surfnodpoint[l]] = m+1;
03632 n++;
03633 }
03634 }
03635 }
03636 }
03637 }
03638
03639
03640
03641
03642
03643 }
03644 }
03645
03646 void fixnodesel::mark_ring_nodes()
03647 {
03648
03649 }
03650
03651
03652
03653
03654
03655
03656
03657
03658
03659
03660
03661
03662
03663
03664
03665
03666 void fixnodesel::select_optimum_2d()
03667 {
03668 double optimum;
03669 double curopt;
03670 long nodes;
03671 long i,j,edges,tnedges;
03672 long stop,e,lnmember,ledges,lnedges;
03673 double zbytek,rat;
03674 bool prime;
03675 bool nondivide;
03676
03677
03678 tnedges = 0;
03679 for(i = 0; i < ncurves; i++){
03680 tnedges+=nedges[i];
03681 }
03682
03683 methodcondcor = curvecond;
03684 typecondcur=n_part_curve;
03685 automember = new long[ncurves];
03686
03687 optimum = 0.05*tnbn;
03688 optimum = round_it(optimum);
03689 if(mespr == 1) fprintf(out,"Optimum is %ld nodes\n",(long)optimum);
03690 for(i = 0; i < ncurves; i++){
03691 nondivide = false;
03692 prime = isPrime(nedges[i]);
03693 switch(nondivide){
03694 case true:
03695 edges=nedges[i]+1;
03696 break;
03697 case false:
03698 edges=nedges[i];
03699 break;
03700 }
03701 rat=(double)nedges[i]/tnedges;
03702 curopt=rat*optimum;
03703 curopt = round_it(curopt);
03704 nodes = (long)curopt;
03705
03706 nodes++;
03707 if(edges%nodes == 0) {
03708 automember[i] = nodes;
03709 }
03710 else{
03711 lnmember = nodes;
03712 ledges = edges;
03713
03714 stop = 0;
03715 e = 0;
03716 while (stop == 0){
03717 zbytek = ledges % lnmember;
03718 if( zbytek != 0){
03719
03720 if(e == 0){
03721 lnmember--;
03722 }
03723 if(e == 1){
03724 lnmember++;
03725 }
03726
03727 stop = 0;
03728 }
03729 else{
03730
03731 stop = 1;
03732 }
03733 if( ledges == lnmember){
03734 e = 1;
03735 lnmember = nmembercur;
03736 stop = 0;
03737 }
03738 if(lnmember == 1){
03739 stop = 1;
03740 e = 2;
03741
03742 ledges++;
03743 lnmember = nmembercur;
03744 }
03745 }
03746 if(e != 2){
03747 automember[i] = lnmember;
03748 }
03749 else{
03750 automember[i] = 1;
03751 }
03752 }
03753 if(mespr == 1) fprintf(out,"curve %ld will be split into %ld parts\n",i,automember[i]);
03754 }
03755
03756
03757
03758
03759
03760
03761
03762
03763
03764
03765
03766
03767
03768
03769
03770
03771
03772
03773
03774
03775
03776
03777
03778
03779
03780
03781
03782
03783
03784
03785
03786
03787
03788
03789
03790
03791
03792
03793
03794
03795
03796
03797
03798
03799
03800
03801
03802
03803
03804
03805
03806
03807
03808
03809
03810
03811
03812
03813
03814
03815
03816
03817
03818
03819
03820
03821
03822
03823
03824
03825
03826
03827
03828
03829
03830
03831
03832
03833
03834
03835
03836
03837
03838
03839
03840
03841
03842
03843
03844
03845
03846
03847
03848
03849 }
03850
03851
03852
03853
03854
03855
03856 void fixnodesel::fixing_detection(long *ltg)
03857 {
03858 double ts,te,tc;
03859 tc = 0.0;
03860 ts = clock();
03861
03862
03863 set_fictitious_subdomain();
03864 te = clock();
03865 tc = tc+te-ts;
03866 if(mespr == 1) fprintf (out,"Time of set of fictitious subdomain %le s\n",(te-ts)/(double)CLOCKS_PER_SEC);
03867
03868 ts = clock();
03869 compute_size_of_fictitious_subdomain();
03870 te = clock();
03871 tc = tc+te-ts;
03872 if(mespr == 1) fprintf (out,"Time of computing of size of fictitious subdomain %le s\n",(te-ts)/(double)CLOCKS_PER_SEC);
03873
03874 ts = clock();
03875 give_whole_dim();
03876 te = clock();
03877 tc = tc+te-ts;
03878 if(mespr == 1) fprintf (out,"Time of set of spatial dimension %le s\n",(te-ts)/(double)CLOCKS_PER_SEC);
03879 fprintf(out,"\n\n\nDimension of problem is %ld\n",dim);
03880 ts = clock();
03881 create_link_dom_nod();
03882 te = clock();
03883 tc = tc+te-ts;
03884 if(mespr == 1) fprintf (out,"Time of creation of globloc link %le s\n",(te-ts)/(double)CLOCKS_PER_SEC);
03885 switch(dim){
03886 case -1:
03887 par_print_err(myrank,"Problem:Spatial dimension is not detected\n", __FILE__, __LINE__, __func__);
03888 break;
03889 case 1:
03890
03891 break;
03892 case 2:
03893
03894 ts = clock();
03895 create_subdom_graph_2d();
03896 te = clock();
03897 tc = tc+te-ts;
03898 if(mespr == 1) fprintf (out,"Time of creation of subdomain graph %le s\n",(te-ts)/(double)CLOCKS_PER_SEC);
03899 ts = clock();
03900 create_master_graph();
03901 te = clock();
03902 tc = tc+te-ts;
03903 if(mespr == 1) fprintf (out,"Time of creation of master graph %le s\n",(te-ts)/(double)CLOCKS_PER_SEC);
03904 ts = clock();
03905 select_minimal_number_2d();
03906 te = clock();
03907 tc = tc+te-ts;
03908 if(mespr == 1) fprintf (out,"Time of selection of minimal number of fixing nodes %le s\n",(te-ts)/(double)CLOCKS_PER_SEC);
03909 if(condfixing == automatic || condfixing == userdef){
03910 if(myrank == 0){
03911 ts = clock();
03912 set_curves();
03913 te = clock();
03914 tc = tc+te-ts;
03915 if(mespr == 1) fprintf (out,"Time of creation of subgraph for adding further fixing nodes %le s\n",(te-ts)/(double)CLOCKS_PER_SEC);
03916 ts = clock();
03917 if(mespr == 1) print_info_curves();
03918 te = clock();
03919 tc = tc+te-ts;
03920 ts = clock();
03921 if(condfixing == automatic){
03922 ts = clock();
03923 select_optimum_2d();
03924 te = clock();
03925 tc = tc+te-ts;
03926 if(mespr == 1) fprintf (out,"Time of selection of optimal number of fixing nodes %le s\n",(te-ts)/(double)CLOCKS_PER_SEC);
03927
03928 }
03929 ts = clock();
03930 add_fixings();
03931 te = clock();
03932 tc = tc+te-ts;
03933 if(mespr == 1) fprintf (out,"Time of adding further fixing nodes %le s\n",(te-ts)/(double)CLOCKS_PER_SEC);
03934 reduce_flag = 0;
03935 }
03936 }
03937 else{
03938 ts = clock();
03939 check_minimal_number_2d();
03940 te = clock();
03941 tc = tc+te-ts;
03942 if(mespr == 1) fprintf (out,"Time of check of minimal number of fixing nodes %le s\n",(te-ts)/(double)CLOCKS_PER_SEC);
03943 }
03944 ts = clock();
03945 assemble_ltg_with_fixings(ltg);
03946 te = clock();
03947 tc = tc+te-ts;
03948 if(mespr == 1) fprintf (out,"Time of assembling of new ltg array %le s\n",(te-ts)/(double)CLOCKS_PER_SEC);
03949 if(mespr == 1) fprintf (out,"Whole time of fixing node detection %le s\n",(tc)/(double)CLOCKS_PER_SEC);
03950 fflush(out);
03951 break;
03952 case 3:
03953
03954 switch(condfixing){
03955 case nocondconer:{
03956 ts = clock();
03957 compute_statistics_of_multiplicity();
03958 te = clock();
03959 tc = tc+te-ts;
03960 if(mespr == 1) fprintf (out,"Time of creation of statistics %le s\n",(te-ts)/(double)CLOCKS_PER_SEC);
03961 minSelStrat = 2;
03962 switch(minSelStrat){
03963 case 1:
03964 if(mespr == 1) fprintf(out,"\nFixing nodes will be select by algorithm based on statictics\n");
03965 if(mespr == 1) fprintf(stdout,"\nFixing nodes will be select by algorithm based on statictics\n");
03966 ts = clock();
03967 select_fixing_nodes_on_master();
03968 te = clock();
03969 tc = tc+te-ts;
03970 if(mespr == 1) fprintf (out,"Time of selection of fixing nodes on master %le s\n",(te-ts)/(double)CLOCKS_PER_SEC);
03971 ts = clock();
03972 check_minimal_number_3d();
03973 te = clock();
03974 tc = tc+te-ts;
03975 if(mespr == 1) fprintf (out,"Time of check of minimal number of fixing nodes %le s\n",(te-ts)/(double)CLOCKS_PER_SEC);
03976 break;
03977 case 2:
03978 if(mespr == 1) fprintf(out,"\nFixing nodes will be select by algortithm based on topology\n");
03979 if(mespr == 1) fprintf(stdout,"\nFixing nodes will be select by algortithm based on topology\n");
03980 ts = clock();
03981 create_subdom_graph_3d ();
03982 te = clock();
03983 tc = tc+te-ts;
03984 if(mespr == 1) fprintf (out,"Time of creation of subdomain graph %le s\n",(te-ts)/(double)CLOCKS_PER_SEC);
03985 ts = clock();
03986 create_master_graph();
03987 te = clock();
03988 tc = tc+te-ts;
03989 if(mespr == 1) fprintf (out,"Time of creation of master graph %le s\n",(te-ts)/(double)CLOCKS_PER_SEC);
03990 ts = clock();
03991 select_minimal_number_3d();
03992 te = clock();
03993 tc = tc+te-ts;
03994 if(mespr == 1) fprintf (out,"Time of selection of minimal number of fixing nodes %le s\n",(te-ts)/(double)CLOCKS_PER_SEC);
03995 ts = clock();
03996 check_minimal_number_3d();
03997 te = clock();
03998 tc = tc+te-ts;
03999 if(mespr == 1) fprintf (out,"Time of check of minimal number of fixing nodes %le s\n",(te-ts)/(double)CLOCKS_PER_SEC);
04000 break;
04001 }
04002 ts = clock();
04003 assemble_ltg_with_fixings(ltg);
04004 te = clock();
04005 tc = tc+te-ts;
04006 if(mespr == 1) fprintf (out,"Time of assembling of new ltg array %le s\n",(te-ts)/(double)CLOCKS_PER_SEC);
04007 if(mespr == 1) fprintf (out,"Whole time of fixing node detection %le s\n",(tc)/(double)CLOCKS_PER_SEC);
04008 fflush(out);
04009 break;
04010 }
04011 case automatic:{
04012
04013 break;
04014 }
04015 case userdef:{
04016 ts = clock();
04017 create_subdom_graph_3d();
04018 te = clock();
04019 tc = tc+te-ts;
04020 if(mespr == 1) fprintf (out,"Time of creation of subdomain graph %le s\n",(te-ts)/(double)CLOCKS_PER_SEC);
04021 ts = clock();
04022 create_master_graph();
04023 te = clock();
04024 tc = tc+te-ts;
04025 if(mespr == 1) fprintf (out,"Time of creation of master graph %le s\n",(te-ts)/(double)CLOCKS_PER_SEC);
04026 ts = clock();
04027 select_minimal_number_3d();
04028 te = clock();
04029 tc = tc+te-ts;
04030 if(mespr == 1) fprintf (out,"Time of selection of minimal number of fixing nodes %le s\n",(te-ts)/(double)CLOCKS_PER_SEC);
04031 ts = clock();
04032 check_minimal_number_3d();
04033 te = clock();
04034 tc = tc+te-ts;
04035 if(mespr == 1) fprintf (out,"Time of check of minimal number of fixing nodes %le s\n",(te-ts)/(double)CLOCKS_PER_SEC);
04036 if(methodcondcor == curvecond || methodcondcor == cursurf) {
04037 if(myrank == 0){
04038 ts = clock();
04039 reduce_master_graph_curve();
04040 te = clock();
04041 tc = tc+te-ts;
04042 if(mespr == 1) fprintf (out,"Time of reduction of master graph to curve graph %le s\n",(te-ts)/(double)CLOCKS_PER_SEC);
04043 ts = clock();
04044 set_curves_3D();
04045 te = clock();
04046 tc = tc+te-ts;
04047 if(mespr == 1) fprintf (out,"Time of creation of subgraph for adding further fixing nodes %le s\n",(te-ts)/(double)CLOCKS_PER_SEC);
04048 ts = clock();
04049 if(mespr == 1) print_info_curves();
04050 te = clock();
04051 tc = tc+te-ts;
04052 ts = clock();
04053 add_fixings();
04054 te = clock();
04055 tc = tc+te-ts;
04056 if(mespr == 1) fprintf (out,"Time of adding further fixing nodes %le s\n",(te-ts)/(double)CLOCKS_PER_SEC);
04057 }
04058 }
04059 if (methodcondcor == surfacecond || methodcondcor == cursurf){
04060 if(myrank == 0){
04061 ts = clock();
04062 reduce_master_graph_surface();
04063 if(mespr == 1) print_info_surfaces();
04064 te = clock();
04065 tc = tc+te-ts;
04066 if(mespr == 1) fprintf (out,"Time of reduction of master graph to surface graph %le s\n",(te-ts)/(double)CLOCKS_PER_SEC);
04067 }
04068 if(typecondsurf == max_tria){
04069 ts = clock();
04070 add_max_triangle_on_surface();
04071 te = clock();
04072 tc = tc+te-ts;
04073 if(mespr == 1) fprintf (out,"Time of adding further fixing nodes into max triangle %le s\n",(te-ts)/(double)CLOCKS_PER_SEC);
04074 }
04075 else{
04076 ts = clock();
04077 select_centers_of_surfaces();
04078 te = clock();
04079 tc = tc+te-ts;
04080 if(mespr == 1) fprintf (out,"Time of selection of surface centres %le s\n",(te-ts)/(double)CLOCKS_PER_SEC);
04081 if(typecondsurf == n_mark || typecondsurf == chose_ring || typecondsurf == chose_max_ring){
04082 if(myrank == 0){
04083 ts = clock();
04084 mark_surfnodes();
04085 te = clock();
04086 tc = tc+te-ts;
04087 }
04088
04089 if(typecondsurf == n_mark || typecondsurf == chose_ring){
04090 order_selected_ring_nodes();
04091 }
04092 }
04093 if(myrank == 0){
04094 ts = clock();
04095 if(mespr == 1) fprintf (out,"Time of selection of marking of surface members %le s\n",(te-ts)/(double)CLOCKS_PER_SEC);
04096 add_fixings();
04097 te = clock();
04098 tc = tc+te-ts;
04099 if(mespr == 1) fprintf (out,"Time of adding further fixing nodes %le s\n",(te-ts)/(double)CLOCKS_PER_SEC);
04100 }
04101 }
04102 }
04103 ts = clock();
04104 assemble_ltg_with_fixings(ltg);
04105 te = clock();
04106 tc = tc+te-ts;
04107 ts = clock();
04108 if(mespr == 1) fprintf (out,"Time of assembling of new ltg array %le s\n",(te-ts)/(double)CLOCKS_PER_SEC);
04109 if(mespr == 1) fprintf (out,"Whole time of fixing node detection %le s\n",(tc)/(double)CLOCKS_PER_SEC);
04110 fflush(out);
04111 break;
04112 }
04113 }
04114 fflush(out);
04115 break;
04116
04117 default:{
04118 par_print_err(myrank,"Spatial dimension is not detected", __FILE__, __LINE__, __func__);
04119 break;
04120 }
04121 }
04122 }
04123
04124
04125
04126
04127
04128
04129
04130
04131
04132
04133
04134 void fixnodesel::reduce_master_graph_curve()
04135 {
04136 long i,j,k,m,l;
04137 long *correduce;
04138
04139 reduce_flag = 3;
04140
04141 ncurnod = 0;
04142 for(i = 0; i < tnbn; i++){
04143 if(multip[i] > 2){
04144 ncurnod++;
04145 }
04146 }
04147
04148 curnod = new long[ncurnod];
04149 curnadjac = new long[ncurnod];
04150 curadjac = new long* [ncurnod];
04151 correduce = new long[tnbn];
04152 curidentif = new long[ncurnod];
04153 k = 0;
04154 l = 0;
04155 for(i = 0; i < tnbn; i++){
04156 if(multip[i] > 2){
04157 curnadjac[k] = 0;
04158 curnod[k] = i;
04159 correduce[i] = k;
04160 curidentif[k] = coarseidentif[i];
04161 for(j = 0; j < coarsenadjac[i]; j++){
04162 if(multip[coarseadjac[i][j]] > 2 ){
04163 curnadjac[k]++;
04164 }
04165 }
04166 curadjac[k] = new long[curnadjac[k]];
04167 m = 0;
04168 for(j = 0; j < coarsenadjac[i]; j++){
04169 if(multip[coarseadjac[i][j]] > 2 ){
04170 curadjac[k][m] = coarseadjac[i][j];
04171 m++;
04172 }
04173 }
04174 k++;
04175 }
04176 else{
04177 correduce[i] = -1;
04178 }
04179 }
04180
04181
04182
04183 for(i = 0; i < ncurnod; i++){
04184 for(j = 0; j < curnadjac[i]; j++){
04185 m = curadjac[i][j];
04186 if(correduce[m] != -1){
04187 curadjac[i][j] = correduce[m];
04188 }
04189 }
04190 }
04191
04192
04193 for( i = 0; i < ncurnod; i++){
04194 if(curnadjac[i] > 2 && curidentif[i] != 3){
04195 curidentif[i] = 3;
04196 }
04197 if(curnadjac[i] == 1 && curidentif[i] != 3){
04198 curidentif[i] = 3;
04199 }
04200 }
04201
04202
04203
04204
04205
04206
04207
04208
04209
04210 delete []correduce;
04211
04212
04213 }
04214
04215
04216
04217
04218
04219
04220
04221
04222 void fixnodesel::reduce_master_graph_surface()
04223 {
04224 long i,j,k,m,l,n;
04225 long a,b;
04226 long auxnsurf;
04227 long *aux;
04228 long **auxsurfmembers;
04229 long *auxnsurfmembers;
04230 long **auxsurfnodmark;
04231
04232 reduce_flag = 2;
04233
04234
04235
04236 subdomsurf = new long*[nproc-1];
04237 for(i = 0; i < nproc-1;i++){
04238 subdomsurf[i] = new long[nproc-i-1];
04239 for(j = 0; j < nproc-i-1; j++){
04240 subdomsurf[i][j] = 0;
04241 }
04242 }
04243
04244 for(i = 0; i < tnbn; i++){
04245 if(multip[i] == 2){
04246 j = glinkdom[i][0];
04247 k = glinkdom[i][1];
04248 if(j < k){
04249 m = k-(j+1);
04250 subdomsurf[j][m]++;
04251 }
04252 else{
04253 m = j-(k+1);
04254 subdomsurf[k][m]++;
04255 }
04256 }
04257 }
04258
04259
04260
04261
04262
04263
04264
04265 auxnsurf = 0;
04266 for(i = 0; i < nproc-1; i++){
04267 for(j = 0; j < nproc-i-1; j++){
04268 if(subdomsurf[i][j] != 0){
04269 auxnsurf++;
04270 }
04271 }
04272 }
04273
04274
04275 auxsurfmembers = new long *[auxnsurf];
04276 auxnsurfmembers = new long [auxnsurf];
04277 m = 0;
04278 for(i = 0; i < nproc-1; i++){
04279 for(j = 0; j < nproc-i-1; j++){
04280 if(subdomsurf[i][j] != 0){
04281 auxsurfmembers[m] = new long[subdomsurf[i][j]];
04282 auxnsurfmembers[m] = 0;
04283 subdomsurf[i][j] = m;
04284 m++;
04285 }
04286 }
04287 }
04288
04289 surfnodpoint = new long[tnbn];
04290 surfnod = new long[tnbn];
04291 n = 0;
04292 for(i = 0; i < tnbn; i++){
04293 if(multip[i] == 2){
04294 j = glinkdom[i][0];
04295 k = glinkdom[i][1];
04296 m = k-(j+1);
04297 l = subdomsurf[j][m];
04298 auxsurfmembers[l][auxnsurfmembers[l]]=i;
04299 surfnodpoint[i]=auxnsurfmembers[l];
04300 surfnod[i]=l;
04301 auxnsurfmembers[l]++;
04302 n++;
04303 }
04304 else{
04305 surfnodpoint[i]= -1;
04306 surfnod[i] = -1;
04307 }
04308 }
04309
04310 for(i = 0; i < nproc-1; i++){
04311 delete []subdomsurf[i];
04312 }
04313 delete []subdomsurf;
04314
04315
04316
04317
04318
04319
04320
04321
04322
04323
04324
04325
04326
04327
04328
04329
04330
04331
04332 l = -2;
04333 aux = new long[auxnsurf];
04334 auxsurfnodmark = new long*[auxnsurf];
04335 for(i = 0; i < auxnsurf; i++){
04336
04337 auxsurfnodmark[i] = new long[auxnsurfmembers[i]];
04338 for(j = 1; j < auxnsurfmembers[i]; j++){
04339 auxsurfnodmark[i][j] = -1;
04340 }
04341 auxsurfnodmark[i][0] = 0;
04342 n = 1;
04343 while(n != auxnsurfmembers[i]){
04344 for(m = 0; m < auxnsurfmembers[i]; m++){
04345 for(j = 0;j < auxnsurfmembers[i]; j++){
04346 if(auxsurfnodmark[i][j] == m){
04347 a = auxsurfmembers[i][j];
04348 for(k = 0; k < coarsenadjac[a]; k++){
04349 b = coarseadjac[a][k];
04350 if(surfnod[b] == i && auxsurfnodmark[i][surfnodpoint[b]] == -1){
04351 auxsurfnodmark[i][surfnodpoint[b]] = m+1;
04352 n++;
04353 }
04354 }
04355 }
04356 if(n == auxnsurfmembers[i]){
04357 break;
04358 }
04359 }
04360 if(n == auxnsurfmembers[i]){
04361 break;
04362 }
04363 }
04364 for(j = 0;j < auxnsurfmembers[i]; j++){
04365 if(auxsurfnodmark[i][j] >= 0){
04366 auxsurfnodmark[i][j] = l;
04367 }
04368 }
04369 aux[i] = -1*l-2;
04370 l--;
04371 if(n != auxnsurfmembers[i]){
04372 for(j = 0;j < auxnsurfmembers[i]; j++){
04373 if(auxsurfnodmark[i][j] == -1){
04374 if(n != auxnsurfmembers[i]-1){
04375 auxsurfnodmark[i][j] = 0;
04376 n++;
04377 }
04378 else{
04379 auxsurfnodmark[i][j] = l;
04380 n++;
04381 aux[i] = -1*l-2;
04382 l--;
04383 }
04384 break;
04385 }
04386 }
04387 }
04388 }
04389
04390
04391
04392
04393
04394 }
04395
04396 nsurf = -1*l-2;
04397
04398 nsurfmembers = new long [nsurf];
04399 surfmembers = new long* [nsurf];
04400 for(i = 0; i < nsurf; i++){
04401 nsurfmembers[i] = 0;
04402 }
04403 if(aux[0] > 1){
04404 for(j = 0; j < auxnsurfmembers[0]; j++){
04405 k = -1*auxsurfnodmark[0][j]-2;
04406 nsurfmembers[k]++;
04407 }
04408 for(j = 0; j < aux[i]+1; j++){
04409 surfmembers[j] = new long[nsurfmembers[j]];
04410 n = 0;
04411 for(k = 0; k < auxnsurfmembers[0]; k++){
04412 a = -1*auxsurfnodmark[0][k]-2;
04413 if(a == j){
04414 b = auxsurfmembers[0][k];
04415 surfnod[b] = j;
04416 surfnodpoint[b] = n;
04417 surfmembers[j][n] = auxsurfmembers[i][k];
04418 n++;
04419 }
04420 }
04421 }
04422 }
04423 else{
04424 nsurfmembers[aux[0]] = auxnsurfmembers[0];
04425 surfmembers[aux[0]] = new long[auxnsurfmembers[0]];
04426 for(k = 0; k < auxnsurfmembers[0]; k++){
04427 b = auxsurfmembers[0][k];
04428 surfnod[b] = aux[0];
04429 surfnodpoint[b] = k;
04430 surfmembers[aux[0]][k] = auxsurfmembers[0][k];
04431 }
04432 }
04433 delete []auxsurfnodmark[0];
04434
04435 for(i = 1; i < auxnsurf; i++){
04436 if(aux[i]-aux[i-1] > 1 ){
04437 for(j = 0; j < auxnsurfmembers[i]; j++){
04438 k = -1*auxsurfnodmark[i][j]-2;
04439 nsurfmembers[k]++;
04440 }
04441 for(j = aux[i-1]+1; j < aux[i]+1; j++){
04442 surfmembers[j] = new long[nsurfmembers[j]];
04443 n = 0;
04444 for(k = 0; k < auxnsurfmembers[i]; k++){
04445 a = -1*auxsurfnodmark[i][k]-2;
04446 if(a == j){
04447 b = auxsurfmembers[i][k];
04448 surfnod[b] = j;
04449 surfnodpoint[b] = n;
04450 surfmembers[j][n] = auxsurfmembers[i][k];
04451 n++;
04452 }
04453 }
04454 }
04455 }
04456 else{
04457 nsurfmembers[aux[i]] = auxnsurfmembers[i];
04458 surfmembers[aux[i]] = new long[nsurfmembers[aux[i]]];
04459 for(k = 0; k < auxnsurfmembers[i]; k++){
04460 b = auxsurfmembers[i][k];
04461 surfnod[b] = aux[i];
04462 surfnodpoint[b] = k;
04463 surfmembers[aux[i]][k] = auxsurfmembers[i][k];
04464 }
04465 }
04466 delete []auxsurfnodmark[i];
04467 }
04468
04469 surfdom = new long*[nsurf];
04470 for(i = 0; i < nsurf; i++){
04471 surfdom[i] = new long[2];
04472 a = surfmembers[i][0];
04473 surfdom[i][0]=glinkdom[a][0];
04474 surfdom[i][1]=glinkdom[a][1];
04475 }
04476
04477
04478 for(i = 0; i < auxnsurf; i++){
04479 delete []auxsurfmembers[i];
04480 }
04481 delete []auxsurfmembers;
04482 delete []auxnsurfmembers;
04483 delete []auxsurfnodmark;
04484
04485
04486 delete []aux;
04487 }
04488
04489
04490
04491
04492
04493
04494
04495
04496 void fixnodesel::set_curves( )
04497 {
04498 long i,j,k;
04499 long cor,comb;
04500 long prev,cur,stop,n,m,a,b,z,v;
04501 long *startnod,*endnod,**corneib,*corpoint;
04502 long **auxmemb;
04503 long *nauxmemb;
04504
04505
04506 cor = 0;
04507 if(coarseidentif != NULL){
04508 for( i = 0; i < tnbn; i++){
04509 if(coarseidentif[i] == 3 ){
04510 cor++;
04511 }
04512 }
04513
04514 if(cor == 0){
04515 for(i = 0; i < tnbn; i++){
04516 if(coarsenadjac[i] == 1 && coarseidentif[i] != 3 ){
04517 coarseidentif[i] = 3;
04518 cor++;
04519 }
04520 if(coarsenadjac[i] > 2 && coarseidentif[i] != 3){
04521 coarseidentif[i] = 3;
04522 cor++;
04523 }
04524 }
04525 }
04526 }
04527 else{
04528 coarseidentif = new long[tnbn];
04529 for(i = 0; i < tnbn; i++){
04530 if(coarsenadjac[i] == 1){
04531 coarseidentif[i] = 3;
04532 }
04533 if(coarsenadjac[i] > 2){
04534 coarseidentif[i] = 3;
04535 }
04536 }
04537 }
04538
04539
04540 corneib = new long*[cor];
04541 corpoint = new long[cor];
04542 k = 0;
04543 comb = 0;
04544 for(i = 0; i < tnbn; i++){
04545 if(coarseidentif[i] == 3 && coarsenadjac[i] != 0){
04546 comb += coarsenadjac[i];
04547 corpoint[k] = i;
04548 corneib[k] = new long[coarsenadjac[i]];
04549 for(j = 0; j < coarsenadjac[i]; j++){
04550 corneib[k][j] = -1;
04551 }
04552 k++;
04553 }
04554 }
04555 if(mespr == 1) fprintf(out,"Estimate number of curves %ld\n",comb);
04556
04557
04558
04559
04560
04561
04562
04563
04564
04565
04566
04567
04568
04569 startnod = new long[comb];
04570 endnod = new long[comb];
04571 auxmemb = new long*[comb];
04572 nauxmemb = new long[comb];
04573
04574
04575 m = 0;
04576 for(i = 0; i < cor; i++){
04577 cur = corpoint[i];
04578 for(j = 0; j < coarsenadjac[corpoint[i]]; j++){
04579 prev = corpoint[i];
04580 if(corneib[i][j] == -1){
04581 startnod[m] = corpoint[i];
04582
04583 stop = 0;
04584 n = 0;
04585 corneib[i][j] = 0;
04586 cur = coarseadjac[corpoint[i]][j];
04587 nauxmemb[m] = 1;
04588 auxmemb[m] = new long[tnbn];
04589 auxmemb[m][n] = cur;
04590 n++;
04591 while(stop != 1){
04592 a = coarsenadjac [cur];
04593 for(k = 0; k < a; k++){
04594 b = coarseadjac[cur][k];
04595 if(b != prev){
04596 if(coarseidentif[b] == 2){
04597 auxmemb[m][n]=b;
04598 n++;
04599 prev = cur;
04600 cur = b;
04601 stop = 0;
04602 }
04603 if(coarseidentif[b] == 3){
04604 stop = 1;
04605 endnod[m] = b;
04606 nauxmemb[m]= n;
04607
04608
04609 for(z = 0; z < cor; z++){
04610 if(b == corpoint[z]){
04611
04612 for(v = 0; v < coarsenadjac[corpoint[z]]; v++){
04613 if(cur == coarseadjac[corpoint[z]][v]){
04614 corneib[z][v] = 0;
04615
04616 break;
04617 }
04618 }
04619 break;
04620 }
04621 }
04622 }
04623 }
04624 }
04625 }
04626 m++;
04627 }
04628 }
04629 }
04630
04631
04632 ncurves = m;
04633 start = new long[ncurves];
04634
04635 end = new long[ncurves];
04636
04637 members = new long*[ncurves];
04638
04639 nmembers = new long[ncurves];
04640
04641 nedges = new long[ncurves];
04642 for(i = 0; i < ncurves; i++){
04643 start[i] = startnod[i];
04644 end[i] = endnod[i];
04645 nmembers[i] = nauxmemb[i];
04646 members[i] = new long[nmembers[i]];
04647 for(j = 0; j < nmembers[i]; j++){
04648 members[i][j] = auxmemb[i][j];
04649 }
04650 nedges[i] = nmembers[i] + 1;
04651 }
04652
04653 delete []startnod;
04654 delete []endnod;
04655 delete []nauxmemb;
04656 for(i = 0; i < ncurves; i++){
04657 delete []auxmemb[i];
04658 }
04659 delete []auxmemb;
04660 delete []corpoint;
04661 for(i = 0; i < cor; i++){
04662 delete []corneib[i];
04663 }
04664 delete []corneib;
04665
04666
04667 }
04668
04669
04670
04671
04672
04673
04674
04675
04676 void fixnodesel::set_curves_3D( )
04677 {
04678 long i,j,k;
04679 long cor,comb;
04680 long prev,cur,stop,n,m,a,b,z,v;
04681 long *startnod,*endnod,**corneib,*corpoint;
04682 long **auxmemb;
04683 long *nauxmemb;
04684
04685 cor = 0;
04686 if(curidentif != NULL){
04687 for( i = 0; i < ncurnod; i++){
04688 if(curidentif[i] == 3 ){
04689 cor++;
04690 }
04691 }
04692
04693 if(cor == 0){
04694 for(i = 0; i < ncurnod; i++){
04695 if(curnadjac[i] == 1){
04696 curidentif[i] = 3;
04697 }
04698 }
04699 }
04700 }
04701 else{
04702 curidentif = new long[ncurnod];
04703 for(i = 0; i < ncurnod; i++){
04704 if(curnadjac[i] == 1){
04705 curidentif[i] = 3;
04706 }
04707 }
04708 }
04709
04710 if(cor == 0){
04711 for(i = 0; i < ncurnod; i++){
04712 if(curnadjac[i] == 1){
04713 curidentif[i] = 3;
04714 }
04715 }
04716 }
04717
04718 cor = 0;
04719 for( i = 0; i < ncurnod; i++){
04720 if(curidentif[i] == 3 ){
04721 cor++;
04722 }
04723 }
04724
04725
04726
04727
04728 corneib = new long*[cor];
04729 corpoint = new long[cor];
04730 k = 0;
04731 comb = 0;
04732 for(i = 0; i < ncurnod; i++){
04733 if(curidentif[i] == 3){
04734 comb += curnadjac[i];
04735 corpoint[k] = i;
04736 corneib[k] = new long[curnadjac[i]];
04737 for(j = 0; j < curnadjac[i]; j++){
04738 corneib[k][j] = -1;
04739 }
04740 k++;
04741 }
04742 }
04743
04744
04745
04746
04747
04748
04749
04750
04751
04752
04753 if(mespr == 1) fprintf(out,"Estimate number of curves %ld\n",comb);
04754
04755 startnod = new long[comb];
04756 endnod = new long[comb];
04757 auxmemb = new long*[comb];
04758 nauxmemb = new long[comb];
04759
04760
04761 m = 0;
04762 for(i = 0; i < cor; i++){
04763 cur = corpoint[i];
04764 for(j = 0; j < curnadjac[corpoint[i]]; j++){
04765 prev = corpoint[i];
04766 if(corneib[i][j] == -1){
04767 startnod[m] = corpoint[i];
04768
04769 stop = 0;
04770 n = 0;
04771 corneib[i][j] = 0;
04772 cur = curadjac[corpoint[i]][j];
04773
04774 if(curidentif[cur] != 3){
04775 nauxmemb[m] = 1;
04776 auxmemb[m] = new long[ncurnod];
04777 auxmemb[m][n] = cur;
04778 n++;
04779 while(stop != 1 ){
04780 a = curnadjac [cur];
04781 for(k = 0; k < a; k++){
04782 b = curadjac[cur][k];
04783 if(b != prev){
04784
04785 if(curidentif[b] == 2){
04786 auxmemb[m][n]=b;
04787 n++;
04788 prev = cur;
04789 cur = b;
04790
04791 stop = 0;
04792
04793 }
04794 if(curidentif[b] == 3){
04795 stop = 1;
04796 endnod[m] = b;
04797 nauxmemb[m]= n;
04798
04799
04800 for(z = 0; z < cor; z++){
04801 if(b == corpoint[z]){
04802
04803 for(v = 0; v < curnadjac[corpoint[z]]; v++){
04804 if(cur == curadjac[corpoint[z]][v]){
04805 corneib[z][v] = 0;
04806
04807 break;
04808 }
04809 }
04810 break;
04811 }
04812 }
04813 break;
04814 }
04815 }
04816 }
04817 }
04818 m++;
04819 }
04820 else{
04821
04822 nauxmemb[m] = 0;
04823 endnod[m] = cur;
04824
04825 for(z = 0; z < cor; z++){
04826 if(cur == corpoint[z]){
04827
04828 for(v = 0; v < curnadjac[corpoint[z]]; v++){
04829 if(prev == curadjac[corpoint[z]][v]){
04830 corneib[z][v] = 0;
04831
04832 break;
04833 }
04834 }
04835 break;
04836 }
04837 }
04838 }
04839 }
04840 }
04841 }
04842
04843
04844 ncurves = m;
04845 start = new long[ncurves];
04846
04847 end = new long[ncurves];
04848
04849 members = new long*[ncurves];
04850
04851 nmembers = new long[ncurves];
04852
04853 nedges = new long[ncurves];
04854 for(i = 0; i < ncurves; i++){
04855 start[i] = startnod[i];
04856 end[i] = endnod[i];
04857 nmembers[i] = nauxmemb[i];
04858 members[i] = new long[nmembers[i]];
04859 for(j = 0; j < nmembers[i]; j++){
04860
04861 k = curnod[auxmemb[i][j]];
04862 members[i][j] = k;
04863 }
04864 nedges[i] = nmembers[i] + 1;
04865 }
04866
04867 delete []startnod;
04868 delete []endnod;
04869 delete []nauxmemb;
04870 for(i = 0; i < ncurves; i++){
04871 delete []auxmemb[i];
04872 }
04873 delete []auxmemb;
04874 delete []corpoint;
04875 for(i = 0; i < cor; i++){
04876 delete []corneib[i];
04877 }
04878 delete []corneib;
04879
04880 }
04881
04882
04883
04884
04885
04886
04887
04888 void fixnodesel::print_info_curves()
04889 {
04890 long i,j;
04891
04892 fprintf(out,"Number of curves %ld\n",ncurves);
04893 for(i = 0; i < ncurves; i++){
04894 fprintf(out,"Curve %ld:\n",i+1);
04895 fprintf(out,"Start %ld\n",start[i]);
04896 fprintf(out,"End %ld\n",end[i]);
04897 fprintf(out,"Number of members %ld\n",nmembers[i]);
04898 fprintf(out,"Members: ");
04899 for(j = 0; j < nmembers[i]; j++){
04900 fprintf(out," %ld",members[i][j]);
04901 }
04902 fprintf(out,"\n");
04903 fprintf(out,"Number of edges %ld\n",nedges[i]);
04904 }
04905 }
04906
04907
04908
04909
04910
04911
04912
04913 void fixnodesel::print_info_surfaces()
04914 {
04915 long i,j;
04916 fprintf(out,"Number of surfaces is %ld\n",nsurf);
04917 for(i = 0; i < nsurf; i++){
04918 fprintf(out,"Surface number %ld lie on subdomain %ld and %ld\n",i,surfdom[i][0],surfdom[i][1]);
04919 fprintf(out,"and has %ld nodes\nMembers: ",nsurfmembers[i]);
04920 for(j = 0; j < nsurfmembers[i]; j++){
04921 fprintf(out," %ld",surfmembers[i][j]);
04922 }
04923 fprintf(out,"\n");
04924 }
04925 }
04926
04927
04928
04929
04930
04931
04932 void fixnodesel::add_fixings()
04933 {
04934
04935 if(methodcondcor == curvecond || methodcondcor == cursurf){
04936 switch(typecondcur){
04937 case nth_memb:{
04938
04939 add_nth_member();
04940 break;
04941 }
04942 case all_memb:{
04943 add_all_members();
04944 break;
04945 }
04946 case rand_memb:{
04947 add_rand();
04948 break;
04949 }
04950 case centroid:{
04951 add_centroid();
04952 break;
04953 }
04954 case n_part_curve:{
04955 add_n_part_curve();
04956 break;
04957 }
04958 case userposdef:{
04959 add_user_pos_def();
04960 break;
04961 }
04962 }
04963 }
04964 if (methodcondcor == surfacecond || methodcondcor == cursurf){
04965
04966 switch(typecondsurf){
04967 case rand_memb:{
04968 add_n_rand_nodes_on_surf();
04969 break;
04970 }
04971 case all_memb:{
04972 add_all_surfnod();
04973 break;
04974 }
04975 case centroid:{
04976 add_centroid_surface();
04977 break;
04978 }
04979 case n_mark:{
04980 add_n_th_mark();
04981 break;
04982 }
04983 case chose_ring:{
04984 add_rings();
04985 break;
04986 }
04987 case chose_max_ring:{
04988 add_max_ring();
04989 break;
04990 }
04991 case userposdef:{
04992 add_user_pos_def();
04993 break;
04994 }
04995 }
04996 }
04997 }
04998
04999
05000
05001
05002
05003
05004 void fixnodesel::add_rand()
05005 {
05006 long i,j,k;
05007 long member;
05008 long *sel_rand;
05009
05010
05011 srand((unsigned int) time(NULL));
05012 sel_rand = new long[nmembercur];
05013 for(i = 0; i < ncurves; i++){
05014 if(mespr == 1) fprintf(out,"On %ld curve will be selected:\n",i);
05015 if(nmembers[i] != 0){
05016 sel_rand[0] = rand()%nmembers[i];
05017 for(j = 1; j < nmembercur; j++){
05018 member = rand()%nmembers[i];
05019
05020
05021
05022
05023
05024
05025
05026
05027
05028
05029
05030
05031
05032
05033
05034
05035 sel_rand[j] = member;
05036 }
05037 for(j = 0; j < nmembercur; j++){
05038 k = members[i][sel_rand[j]];
05039 coarseidentif[k] = 3;
05040 if(mespr == 1) fprintf(out," %ld",members[i][sel_rand[j]]);
05041 }
05042 }
05043 if(mespr == 1) fprintf(out,"\n");
05044 }
05045
05046
05047 delete []sel_rand;
05048 fflush(out);
05049
05050 }
05051
05052
05053
05054
05055
05056
05057 void fixnodesel::add_each_nth_member()
05058 {
05059 long i,j,stop,e;
05060 double zbytek;
05061 long cycles,member,lnmember;
05062
05063 for(i = 0; i < ncurves; i++){
05064 if(mespr == 1) fprintf(out,"On curve %ld will be selected:",i);
05065 lnmember = nmembercur;
05066 if(nmembers[i] != 0){
05067 if(lnmember > nedges[i]){
05068
05069 lnmember = nedges[i]-1;
05070 }
05071 zbytek = nedges[i] % lnmember;
05072 if(zbytek == 0){
05073 cycles = nedges[i] / lnmember;
05074 }
05075 else{
05076
05077 stop = 0;
05078 e = 0;
05079 while (stop == 0){
05080 zbytek = nedges[i] % lnmember;
05081 if( zbytek != 0){
05082
05083 if(e == 0){
05084 lnmember++;
05085 }
05086 if(e == 1){
05087 lnmember--;
05088 }
05089
05090 stop = 0;
05091 }
05092 else{
05093 cycles = nedges[i] / lnmember;
05094
05095 stop = 1;
05096 }
05097 if( nedges[i] == lnmember){
05098 e = 1;
05099 lnmember = nmembercur;
05100 stop = 0;
05101 }
05102 if(lnmember == 1){
05103 stop = 1;
05104 e = 2;
05105 }
05106 }
05107 }
05108 if(e == 2){
05109
05110 member = 0;
05111 for(j = 0; j < nmembers[i]; j++){
05112 member++;
05113 if(member == nmembercur){
05114
05115 member = 0;
05116 coarseidentif[members[i][j]] = 3;
05117 }
05118 }
05119 }
05120
05121 for(j = 1; j < cycles; j++){
05122 member = j*lnmember-1;
05123 if(mespr == 1) fprintf(out," %ld",members[i][member]);
05124 coarseidentif[members[i][member]] = 3;
05125 }
05126 }
05127 if(mespr == 1) fprintf(out,"\n");
05128 }
05129 fflush(out);
05130 }
05131
05132
05133
05134
05135
05136
05137 void fixnodesel::add_nth_member()
05138 {
05139 long i,j;
05140 long member;
05141
05142 for(i = 0; i < ncurves; i++){
05143 if(condfixing == automatic){
05144 nmembercur = automember[i];
05145 }
05146
05147 if(mespr == 1) fprintf(out,"On curve %ld will be selected:",i);
05148 if(nmembers[i] != 0){
05149 if(nmembercur >= nmembers[i]){
05150
05151 }
05152 else{
05153 member = 0;
05154 for(j = 0; j < nmembers[i]; j++){
05155 member++;
05156 if(member == nmembercur){
05157 if(mespr == 1) fprintf(out," %ld",members[i][j]);
05158 member = 0;
05159 coarseidentif[members[i][j]] = 3;
05160 }
05161 }
05162 }
05163 }
05164 if(mespr == 1) fprintf(out,"\n");
05165 }
05166 fflush(out);
05167 }
05168
05169
05170
05171
05172
05173
05174
05175 void fixnodesel::add_all_members( )
05176 {
05177 long i,j;
05178 for(i = 0; i < ncurves; i++){
05179 if(mespr == 1) fprintf(out,"On curve %ld will be selected:",i);
05180 for(j = 0; j < nmembers[i]; j++){
05181 if(mespr == 1) fprintf(out," %ld",members[i][j]);
05182 coarseidentif[members[i][j]] = 3;
05183 }
05184 if(mespr == 1) fprintf(out,"\n");
05185 }
05186 }
05187
05188
05189
05190
05191
05192
05193
05194 void fixnodesel:: add_centroid()
05195 {
05196 long i,k;
05197 long select;
05198 for(i = 0; i < ncurves; i++){
05199 if((nedges[i]%2) == 0){
05200 select = nedges[i]/2;
05201 select--;
05202 }
05203 else{
05204 k = nedges[i]+1;
05205 select = k/2;
05206 select--;
05207 }
05208 if(mespr == 1) fprintf(out,"On curve %ld will be selected %ld\n",i,members[i][select]);
05209 coarseidentif[members[i][select]] = 3;
05210 }
05211
05212 }
05213
05214
05215
05216
05217
05218
05219
05220
05221
05222
05223
05224
05225
05226
05227
05228
05229
05230
05231
05232
05233
05234
05235
05236
05237
05238
05239
05240
05241
05242
05243
05244
05245
05246
05247
05248
05249
05250
05251
05252
05253
05254
05255
05256
05257
05258
05259
05260
05261
05262
05263
05264
05265
05266
05267
05268
05269
05270
05271
05272
05273
05274
05275
05276
05277
05278
05279
05280
05281
05282
05283
05284
05285
05286
05287
05288
05289
05290
05291
05292
05293
05294
05295
05296
05297
05298
05299
05300
05301
05302
05303
05304
05305
05306
05307
05308
05309
05310
05311
05312
05313
05314
05315
05316
05317
05318
05319
05320
05321
05322
05323
05324
05325
05326
05327
05328
05329
05330
05331
05332
05333
05334
05335
05336
05337
05338
05339
05340
05341
05342
05343
05344
05345
05346
05347
05348
05349
05350
05351
05352
05353
05354
05355
05356
05357
05358
05359
05360
05361
05362
05363
05364
05365
05366
05367
05368
05369
05370
05371
05372
05373
05374
05375
05376
05377
05378
05379
05380
05381
05382
05383
05384
05385
05386
05387
05388
05389
05390
05391
05392
05393
05394
05395
05396
05397
05398
05399
05400
05401
05402
05403
05404
05405
05406
05407
05408
05409
05410
05411
05412
05413
05414
05415
05416
05417
05418
05419
05420
05421
05422
05423
05424
05425
05426
05427
05428
05429
05430
05431
05432
05433
05434
05435
05436
05437
05438
05439
05440
05441
05442
05443
05444
05445
05446
05447
05448
05449
05450
05451
05452
05453
05454
05455
05456
05457
05458
05459
05460
05461
05462
05463
05464
05465
05466
05467
05468
05469
05470
05471
05472
05473
05474
05475
05476
05477
05478
05479
05480
05481
05482
05483
05484
05485
05486
05487
05488
05489
05490
05491
05492
05493
05494
05495
05496
05497
05498
05499
05500
05501
05502
05503
05504
05505
05506
05507
05508
05509
05510
05511
05512
05513
05514
05515
05516
05517
05518
05519
05520
05521
05522
05523
05524
05525
05526
05527
05528
05529
05530
05531
05532
05533
05534
05535
05536
05537
05538
05539
05540
05541
05542 void fixnodesel::add_n_part_curve()
05543 {
05544 long i,j,k;
05545 long smember,stop,e,lnmember,lnedges;
05546 double zbytek;
05547
05548 for(i = 0; i < ncurves; i++){
05549 if(condfixing == automatic){
05550 nmembercur = automember[i];
05551 }
05552 if(mespr == 1) fprintf(out,"On curve %ld will be selected:",i);
05553 if(nedges[i] > nmembercur ){
05554 if(nedges[i]%nmembercur == 0) {
05555 smember = nedges[i]/nmembercur;
05556 }
05557 else{
05558 lnmember = nmembercur;
05559 lnedges = nedges[i];
05560
05561 stop = 0;
05562 e = 0;
05563 while (stop == 0){
05564 zbytek = lnedges % lnmember;
05565 if( zbytek != 0){
05566
05567 if(e == 0){
05568 lnmember--;
05569 }
05570 if(e == 1){
05571 lnmember++;
05572 }
05573
05574 stop = 0;
05575 }
05576 else{
05577 smember = lnedges / lnmember;
05578
05579 stop = 1;
05580 }
05581 if( lnedges == lnmember){
05582 e = 1;
05583 lnmember = nmembercur;
05584 stop = 0;
05585 }
05586 if(lnmember == 1){
05587 stop = 1;
05588 e = 2;
05589
05590 lnedges++;
05591 lnmember = nmembercur;
05592 }
05593 }
05594 }
05595
05596 if(e != 2){
05597 k = 0;
05598 for(j = 0; j < nmembers[i];j++){
05599 k++;
05600 if(k == smember){
05601 k = 0;
05602 if(mespr == 1) fprintf(out," %ld",members[i][j]);
05603 coarseidentif[members[i][j]] = 3;
05604 }
05605 }
05606 }
05607 }
05608 if(mespr == 1) fprintf(out,"\n");
05609 }
05610
05611 }
05612
05613
05614
05615
05616
05617 void fixnodesel::add_user_pos_def()
05618 {
05619 long i,j,k,select,numbering;
05620 if(methodcondcor == curvecond || methodcondcor == cursurf){
05621 numbering =1;
05622 for(i = 0; i < ncurves; i++){
05623 if(mespr == 1) fprintf(out,"On curve %ld will be selected\n",i);
05624 coarseidentif[start[i]] = 2;
05625 coarseidentif[end[i]] = 2;
05626 for(j = 0; j < nmembers[i]; j++){
05627 coarseidentif[members[i][j]] = 2;
05628 }
05629
05630 for(j = 0; j < nuserdefnod; j++){
05631 select = userdefnod[j]-1;
05632 if(select == -1 || numbering == 0){
05633 select++;
05634 numbering = 0;
05635 }
05636
05637 if(select < nmembers[i]){
05638 k = members[i][select];
05639 coarseidentif[k] = 3;
05640 if(mespr == 1) fprintf(out," %ld",k);
05641 }
05642 else{
05643 select = nmembers[i]-1;
05644 k = members[i][select];
05645 coarseidentif[k] = 3;
05646 if(mespr == 1) fprintf(out," %ld",k);
05647 break;
05648 }
05649 }
05650 if(mespr == 1) fprintf(out,"\n");
05651 }
05652 }
05653 if (methodcondcor == surfacecond || methodcondcor == cursurf){
05654 numbering =1;
05655 for(i = 0; i < nsurf; i++){
05656 if(mespr == 1) fprintf(out,"On surface %ld will be selected\n",i);
05657 for(j = 0; j < nsurfmembers[i]; j++){
05658 coarseidentif[surfmembers[i][j]] = 2;
05659 }
05660
05661 for(j = 0; j < nuserdefnod; j++){
05662 select = userdefnod[j]-1;
05663 if(select == -1 || numbering == 0){
05664 select++;
05665 numbering = 0;
05666 }
05667
05668 if(select < nsurfmembers[i]){
05669 k = surfmembers[i][select];
05670 coarseidentif[k] = 3;
05671 if(mespr == 1) fprintf(out," %ld",k);
05672 }
05673 else{
05674 select = nmembers[i]-1;
05675 k = surfmembers[i][select];
05676 coarseidentif[k] = 3;
05677 if(mespr == 1) fprintf(out," %ld",k);
05678 break;
05679 }
05680 }
05681 if(mespr == 1) fprintf(out,"\n");
05682 }
05683 }
05684
05685 }
05686
05687
05688
05689
05690
05691 void fixnodesel::add_n_rand_nodes_on_surf()
05692 {
05693
05694 long i,j,k;
05695 long member;
05696 long *sel_rand;
05697
05698
05699 srand((unsigned int) time(NULL));
05700 sel_rand = new long[nmembersurf];
05701 for(i = 0; i < nsurf; i++){
05702 if(mespr == 1) fprintf(out,"On %ld surface will be selected:\n",i);
05703 if(nsurfmembers[i] != 0){
05704 sel_rand[0] = rand()%nsurfmembers[i];
05705 for(j = 1; j < nmembersurf; j++){
05706 member = rand()%nsurfmembers[i];
05707
05708
05709
05710
05711
05712
05713
05714
05715
05716
05717
05718
05719
05720
05721
05722
05723 sel_rand[j] = member;
05724 }
05725 for(j = 0; j < nmembersurf; j++){
05726 k = surfmembers[i][sel_rand[j]];
05727 coarseidentif[k] = 3;
05728 if(mespr == 1) fprintf(out," %ld",k);
05729 }
05730 }
05731 if(mespr == 1) fprintf(out,"\n");
05732 }
05733
05734
05735 delete []sel_rand;
05736 fflush(out);
05737 }
05738
05739
05740
05741
05742
05743 void fixnodesel::add_all_surfnod()
05744 {
05745 long i,j;
05746 for(i = 0; i < nsurf; i++){
05747 if(mespr == 1) fprintf(out,"On surface %ld will be selected:",i);
05748 for(j = 0; j < nsurfmembers[i]; j++){
05749 if(mespr == 1) fprintf(out," %ld",surfmembers[i][j]);
05750 coarseidentif[surfmembers[i][j]] = 3;
05751 }
05752 if(mespr == 1) fprintf(out,"\n");
05753 }
05754 }
05755
05756
05757
05758
05759
05760
05761
05762 void fixnodesel::add_centroid_surface()
05763 {
05764 long i,j,a;
05765 for(i = 0; i < nsurf; i++){
05766 if(mespr == 1) fprintf(out,"On surface %ld will be selected:",i);
05767 for(j = 0; j < nsurfcentres[i]; j++){
05768 a = surfcenters[i][j];
05769 coarseidentif[a] = 3;
05770 if(mespr == 1) fprintf(out," %ld",a);
05771 }
05772 if(mespr == 1) fprintf(out,"\n");
05773 }
05774
05775 }
05776
05777
05778
05779
05780
05781
05782 void fixnodesel::add_n_th_mark()
05783 {
05784 long i,j,a;
05785 for(i = 0; i < nsurf; i++){
05786 if(mespr == 1) fprintf(out,"On surface %ld will be selected:",i);
05787 if(nsurfmembers[i] != 0){
05788 for(j = 0; j < nsurfmembers[i]; j++){
05789 if((surfnodmark[i][j]%nmembersurf) == 0){
05790 a = surfmembers[i][j];
05791 coarseidentif[a] = 3;
05792 if(mespr == 1) fprintf(out," %ld",surfmembers[i][j]);
05793 }
05794 }
05795 }
05796 if(mespr == 1) fprintf(out,"\n");
05797 }
05798 }
05799
05800
05801
05802
05803
05804 void fixnodesel::add_int_points_surface()
05805 {
05806
05807 }
05808
05809
05810
05811
05812
05813 void fixnodesel::order_selected_ring_nodes()
05814 {
05815 long i,j,a,k,max,m;
05816 long nring;
05817 long *ringnodes;
05818 long *ringmark;
05819 long **pointer;
05820 bool stop;
05821 MPI_Status stat;
05822
05823
05824 ringnodes = new long[nbn];
05825
05826 if(myrank == 0){
05827 pointer = new long*[nproc];
05828 for(i = 0; i< nproc; i++){
05829 pointer[i] = new long[nbn];
05830 for (j = 0; j < nbn; j++){
05831 pointer[i][j] = -1;
05832 }
05833 }
05834 long aa,bb,n;
05835 for(i = 0; i < nsurf; i++){
05836 max = 0;
05837 aa = surfdom[i][0];
05838 bb = surfdom[i][1];
05839 for(j = 0; j < nsurfmembers[i]; j++){
05840 if(max < surfnodmark[i][j]){
05841 max = surfnodmark[i][j];
05842 }
05843 }
05844 if((max%nmembersurf) == 0){
05845 nring = max/nmembersurf;
05846 }
05847 a = nmembersurf;
05848 for (j = 0; j < nring; j++){
05849 for(k = 0; k < nsurfmembers[i]; k++){
05850 if(surfnodmark[i][k] == a){
05851
05852 pointer[aa][surfnodpoint[surfmembers[i][k]]] = i;
05853 pointer[bb][surfnodpoint[surfmembers[i][k]]] = i;
05854 }
05855 a+=nmembersurf;
05856 }
05857 }
05858 }
05859
05860 for(i = 1; i < nproc; i++){
05861 MPI_Send (pointer[i],nbnd[i],MPI_LONG,i,myrank,MPI_COMM_WORLD);
05862 }
05863 for(i = 0; i < nbnd[0]; i++){
05864 ringnodes[i] = pointer[0][i];
05865 }
05866 }
05867 else{
05868
05869 MPI_Recv (ringnodes,nbn,MPI_LONG,MPI_ANY_SOURCE,MPI_ANY_TAG,MPI_COMM_WORLD,&stat);
05870 }
05871 MPI_Barrier (MPI_COMM_WORLD);
05872
05873
05874 long *aux;
05875 max = 0;
05876 for(i = 0; i < nbn; i++){
05877 if(max < ringnodes[i]){
05878 max = ringnodes[i];
05879 }
05880 }
05881
05882 max++;
05883 aux = new long[max];
05884 for(i = 0; i < max; i++){
05885 aux[i] = 0;
05886 }
05887
05888 for(i = 0; i < nbn; i++){
05889 if(ringnodes[i] != -1){
05890 j = ringnodes[i];
05891 aux[j]++;
05892 }
05893 }
05894 double *nodelength;
05895 double *nodeangle;
05896 nodelength = new double[nbn];
05897 nodeangle = new double[nbn];
05898 k = 0;
05899 for(i = 0; i < max; i++){
05900 if(aux[i] != 0){
05901 for(j = 0; j < nbn; j++){
05902 if(ringnodes[j] == i){
05903 nodelength[j] = (top->gnodes[lnbn[j]].x-realcg[k][0])*(top->gnodes[lnbn[j]].x-realcg[k][0])+(top->gnodes[lnbn[j]].y-realcg[k][1])*(top->gnodes[lnbn[j]].y-realcg[k][1])+(top->gnodes[lnbn[j]].z-realcg[k][2])*(top->gnodes[lnbn[j]].z-realcg[k][2]);
05904 nodeangle[j] = (top->gnodes[lnbn[j]].x-realcg[k][0])*(top->gnodes[lnbn[j]].x-realcg[k][0])+(top->gnodes[lnbn[j]].y-realcg[k][1])*(top->gnodes[lnbn[j]].y-realcg[k][1])+(top->gnodes[lnbn[j]].z-realcg[k][2])*(top->gnodes[lnbn[j]].z-realcg[k][2]);
05905 }
05906 }
05907
05908
05909 k++;
05910 }
05911 }
05912
05913
05914 }
05915
05916
05917
05918
05919
05920
05921 void fixnodesel::add_rings()
05922 {
05923 long i,j,k,a;
05924 for(i = 0; i < nsurf; i++){
05925 if(mespr == 1) fprintf(out,"On surface %ld will be selected:",i);
05926 if(nring != 0){
05927 for(j = 0; j < nring; j++){
05928 for(k = 0; k < nsurfmembers[i]; k++){
05929 if(surfnodmark[i][k] == ring[j]){
05930 a = surfmembers[i][j];
05931 coarseidentif[a] = 3;
05932 if(mespr == 1) fprintf(out," %ld",surfmembers[i][k]);
05933 }
05934 }
05935 if(mespr == 1) fprintf(out,"\n");
05936 }
05937 }
05938 if(mespr == 1) fprintf(out,"\n");
05939 }
05940
05941 }
05942
05943
05944
05945
05946 void fixnodesel::add_max_ring()
05947 {
05948 long i,j,a;
05949 long max;
05950 for(i = 0; i < nsurf; i++){
05951 max = 0;
05952 for(j = 0; j < nsurfmembers[i]; j++){
05953 if(max < surfnodmark[i][j]){
05954 max = surfnodmark[i][j];
05955 }
05956 }
05957 if(mespr == 1) fprintf(out,"On surface %ld will be selected:",i);
05958 for(j = 0; j < nsurfmembers[i]; j++){
05959 if(surfnodmark[i][j] == max || surfnodmark[i][j] == max-1){
05960 a = surfmembers[i][j];
05961 coarseidentif[a] = 3;
05962 if(mespr == 1) fprintf(out," %ld",surfmembers[i][j]);
05963 }
05964 }
05965 if(mespr == 1) fprintf(out,"\n");
05966 }
05967
05968 }
05969
05970
05971
05972
05973
05974
05975 void fixnodesel::add_max_triangle_on_surface()
05976 {
05977 long i,j,k,l,max;
05978 long buffsize,*buff,*aux;
05979 MPI_Status stat;
05980
05981
05982 double maxradius,radius,maxarea,area;
05983 long **auxnodes;
05984 long *vcomb;
05985 long subdomnsurf;
05986 long *subdomnsurfmem;
05987 long **subdomsurfmem;
05988
05989 buffsize = maxnbn;
05990 buff = new long[buffsize];
05991 if(myrank == 0){
05992 for(i = 1; i < nproc; i++){
05993 for(j = 0; j < nbnd[i]; j++){
05994 buff[j] = surfnod[cnbn[i][j]];
05995 }
05996 MPI_Send (buff,buffsize,MPI_LONG,i,myrank,MPI_COMM_WORLD);
05997
05998 }
05999 for(j = 0; j < nbnd[0]; j++){
06000 buff[j] = surfnod[cnbn[0][j]];
06001 }
06002 }
06003 else{
06004 MPI_Recv (buff,buffsize,MPI_LONG,MPI_ANY_SOURCE,MPI_ANY_TAG,MPI_COMM_WORLD,&stat);
06005 }
06006 MPI_Barrier (MPI_COMM_WORLD);
06007
06008 max = 0;
06009 for(i = 0; i < nbn; i++){
06010 if(max < buff[i]){
06011 max = buff[i];
06012 }
06013 }
06014
06015
06016 max++;
06017 aux = new long[max];
06018 for(i = 0; i < max; i++){
06019 aux[i] = 0;
06020 }
06021
06022 for(i = 0; i < nbn; i++){
06023 if(buff[i] != -1){
06024 j = buff[i];
06025 aux[j]++;
06026 }
06027 }
06028
06029 subdomnsurf = 0;
06030 for(i = 0; i < max; i++){
06031 if(aux[i] != 0){
06032 subdomnsurf++;
06033 }
06034 }
06035
06036 subdomnsurfmem = new long[subdomnsurf];
06037 subdomsurfmem = new long*[subdomnsurf];
06038 for(i = 0; i < subdomnsurf; i++){
06039 subdomnsurfmem[i] = 0;
06040 }
06041
06042 j = 0;
06043 for(i = 0; i < max; i++){
06044 if(aux[i] != 0){
06045 subdomsurfmem[j] = new long[aux[i]];
06046 aux[i] = j;
06047 j++;
06048 }
06049 }
06050
06051 for(i = 0; i < nbn; i++){
06052 if(buff[i] != -1){
06053 j = buff[i];
06054 k = aux[j];
06055 subdomsurfmem[k][subdomnsurfmem[k]]=i;
06056 subdomnsurfmem[k]++;
06057 buff[i] = 0;
06058 }
06059 }
06060 delete[]aux;
06061
06062
06063
06064
06065
06066
06067
06068
06069
06070
06071
06072
06073 max = 0;
06074 for(i = 0; i < subdomnsurf; i++){
06075 if(max < subdomnsurfmem[i]){
06076 max = subdomnsurfmem[i];
06077 }
06078 }
06079
06080 vcomb = new long[3];
06081 for(i = 0; i < subdomnsurf; i++){
06082 vcomb[0] = lnbn[subdomsurfmem[i][0]];
06083 maxradius = 0.0;
06084 for(j = 0; j < subdomnsurfmem[i]; j++){
06085 vcomb[1] = lnbn[subdomsurfmem[i][j]];
06086 radius = compute_length_of_vector(vcomb[0],vcomb[1]);
06087
06088 if(radius > maxradius){
06089 maxradius = radius;
06090 k = j;
06091 }
06092 }
06093 vcomb[0] = lnbn[subdomsurfmem[i][k]];
06094
06095
06096 maxradius = 0.0;
06097 for(j = 0; j < subdomnsurfmem[i]; j++){
06098 vcomb[1] = lnbn[subdomsurfmem[i][j]];
06099 radius = compute_length_of_vector(vcomb[0],vcomb[1]);
06100 if(radius > maxradius){
06101 k = j;
06102 maxradius = radius;
06103 }
06104 }
06105 vcomb[0] = lnbn[subdomsurfmem[i][k]];
06106 buff[subdomsurfmem[i][k]]=1;
06107
06108 maxradius = 0.0;
06109 for(j = 0; j < subdomnsurfmem[i]; j++){
06110 vcomb[1] = lnbn[subdomsurfmem[i][j]];
06111 radius = compute_length_of_vector(vcomb[0],vcomb[1]);
06112 if(radius > maxradius){
06113 k = j;
06114 maxradius = radius;
06115 }
06116 }
06117 vcomb[1] = lnbn[subdomsurfmem[i][k]];
06118 buff[subdomsurfmem[i][k]]=1;
06119
06120
06121 maxarea = 0.0;
06122 for(j = 0; j < subdomnsurfmem[i]; j++){
06123 if(vcomb[0] != lnbn[subdomsurfmem[i][j]] && vcomb[1] != lnbn[subdomsurfmem[i][j]]){
06124 vcomb[2] = lnbn[subdomsurfmem[i][j]];
06125 area = compute_area_of_triangle(vcomb[0],vcomb[1],vcomb[2]);
06126 if(maxarea < area){
06127 maxarea = area;
06128 k = j;
06129 }
06130 }
06131 }
06132 buff[subdomsurfmem[i][k]]=1;
06133
06134 }
06135 delete []vcomb;
06136
06137
06138 if(myrank == 0){
06139 auxnodes = new long*[nproc];
06140 auxnodes[0] = new long[nbnd[0]];
06141 for(j = 0; j < nbnd[0]; j++){
06142 auxnodes[0][j] = buff[j];
06143 }
06144 for(i = 1; i < nproc; i++){
06145 MPI_Recv (buff,buffsize,MPI_LONG,MPI_ANY_SOURCE,MPI_ANY_TAG,MPI_COMM_WORLD,&stat);
06146 k = domproc[stat.MPI_TAG];
06147 auxnodes[k] = new long[nbnd[k]];
06148 for(j = 0; j < nbnd[k]; j++){
06149 auxnodes[k][j] = buff[j];
06150 }
06151 }
06152 }
06153 else{
06154 MPI_Send (buff,buffsize,MPI_LONG,0,myrank,MPI_COMM_WORLD);
06155 }
06156 MPI_Barrier (MPI_COMM_WORLD);
06157 delete []buff;
06158
06159 if(myrank == 0){
06160 for(i = 0; i < nproc; i++){
06161 for(j = 0; j < nbnd[i]; j++){
06162 if(auxnodes[i][j] == 1){
06163 k = cnbn[i][j];
06164 l = surfnod[k];
06165 if(glinkdom[k][1] <= i){
06166 coarseidentif[k] = 3;
06167 }
06168 }
06169 }
06170 }
06171
06172 for(i = 0; i < nsurf; i++){
06173 if(mespr == 1) fprintf(out,"On surface %ld will be selected:",i);
06174 for(j = 0; j < nsurfmembers[i]; j++){
06175 k = surfmembers[i][j];
06176 if(coarseidentif[k] == 3){
06177 if(mespr == 1) fprintf(out," %ld",surfmembers[i][j]);
06178 }
06179 }
06180 if(mespr == 1) fprintf(out,"\n");
06181 }
06182 for(i = 0; i < nproc; i++){
06183 delete []auxnodes[i];
06184 }
06185 delete []auxnodes;
06186 }
06187 }
06188
06189
06190