00001 #include "partopjk.h"
00002 #include "psolver.h"
00003 #include <stdlib.h>
00004 #include <string.h>
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017 partopjk::partopjk(int np,int mr,meshdescription meshd,char *nameproc, int nl)
00018 {
00019
00020 nproc=np;
00021
00022 myrank=mr;
00023
00024 md = meshd;
00025
00026 strcpy (procName,nameproc);
00027 nameLength=nl;
00028
00029
00030 nn = 0;
00031
00032 maxnn = 0;
00033
00034 tnnp = 0;
00035
00036 maxndofd = 0;
00037
00038 nin = 0;
00039
00040 nbn = 0;
00041
00042 nidof = 0;
00043
00044 nbdof = 0;
00045
00046 ndof = 0;
00047
00048 maxnbdofd = 0;
00049
00050
00051
00052
00053
00054
00055 tnnp = 0;
00056
00057 ndofc = 0;
00058
00059
00060
00061
00062
00063
00064 bnid=NULL;
00065
00066 inid=NULL;
00067
00068
00069
00070
00071
00072
00073 nnsd = NULL;
00074
00075 mltg = NULL;
00076
00077 nodmultip = NULL;
00078
00079 noddom = NULL;
00080
00081 ndofnm = NULL;
00082
00083 dofm = NULL;
00084
00085 nbnd = NULL;
00086
00087 nind = NULL;
00088
00089 nbdofd = NULL;
00090
00091 bdofd = NULL;
00092
00093
00094 fullnoddom = NULL;
00095
00096 lagrmultip = NULL;
00097 }
00098
00099
00100
00101
00102
00103
00104 partopjk::~partopjk()
00105 {
00106 long i;
00107
00108 if (myrank==0){
00109
00110 delete [] nnsd;
00111
00112
00113 for (i=0;i<nproc;i++){
00114 delete [] mltg[i];
00115 }
00116 delete [] mltg;
00117
00118
00119 delete [] nodmultip;
00120
00121
00122 delete [] noddom;
00123
00124
00125 delete [] ndofnm;
00126
00127
00128 for (i=0;i<tnnp;i++){
00129 delete [] dofm[i];
00130 }
00131 delete [] dofm;
00132
00133
00134 delete [] nbnd;
00135
00136 delete [] nind;
00137
00138 delete [] nbdofd;
00139
00140
00141 for (i=0;i<nproc;i++){
00142 delete [] bdofd[i];
00143 }
00144 delete [] bdofd;
00145
00146
00147
00148 for (i=0;i<tnnp;i++){
00149 delete [] fullnoddom[i];
00150 }
00151 delete [] fullnoddom;
00152
00153
00154 for (i=0;i<tnnp;i++){
00155 delete [] lagrmultip[i];
00156 }
00157 delete [] lagrmultip;
00158
00159 }
00160
00161
00162 delete [] bnid;
00163
00164 delete [] inid;
00165 }
00166
00167
00168
00169
00170
00171
00172
00173
00174
00175
00176
00177 void partopjk::numbers_of_nodes_on_subdomains (gtopology *top,long *domproc,FILE *out)
00178 {
00179 long i,j,k;
00180 MPI_Status stat;
00181
00182
00183 nn = top->nn;
00184
00185 if (myrank==0){
00186 if (nnsd!=NULL){
00187 delete [] nnsd;
00188 }
00189 nnsd = new long [nproc];
00190
00191
00192 j=domproc[0];
00193 nnsd[j]=nn;
00194 if (maxnn<nn){
00195 maxnn=nn;
00196 }
00197
00198
00199 for (i=1;i<nproc;i++){
00200 MPI_Recv (&k,1,MPI_LONG,MPI_ANY_SOURCE,MPI_ANY_TAG,MPI_COMM_WORLD,&stat);
00201 j=domproc[stat.MPI_TAG];
00202 nnsd[j]=k;
00203 if (maxnn<k){
00204 maxnn=k;
00205 }
00206 }
00207
00208 for (i=1;i<nproc;i++){
00209 MPI_Send (&maxnn,1,MPI_LONG,i,myrank,MPI_COMM_WORLD);
00210 }
00211
00212 }
00213 else{
00214 MPI_Send (&nn,1,MPI_LONG,0,myrank,MPI_COMM_WORLD);
00215 MPI_Recv (&maxnn,1,MPI_LONG,MPI_ANY_SOURCE,MPI_ANY_TAG,MPI_COMM_WORLD,&stat);
00216 }
00217 MPI_Barrier (MPI_COMM_WORLD);
00218
00219
00220 fprintf (out,"\n\n maxnn (maximum number of nodes on subdomain) %ld",maxnn);
00221 if (myrank==0){
00222 fprintf (out,"\n\n array nnsd (number of nodes on subdomains");
00223 for (i=0;i<nproc;i++){
00224 fprintf (out,"\n nnsd %6ld %ld",i+1,nnsd[i]);
00225 }
00226 }
00227
00228 }
00229
00230
00231
00232
00233
00234
00235
00236
00237
00238
00239 void partopjk::ltg_on_master (long *ltg,long *domproc,FILE *out)
00240 {
00241 long i,j,k,buffsize;
00242 long *buff;
00243 MPI_Status stat;
00244
00245 buffsize=maxnn;
00246 buff = new long [buffsize];
00247
00248 for (i=0;i<nn;i++){
00249 buff[i]=ltg[i];
00250 }
00251
00252 if (myrank==0){
00253
00254 mltg = new long* [nproc];
00255 for (i=0;i<nproc;i++){
00256 mltg[i] = new long [nnsd[i]];
00257 }
00258
00259
00260 j=domproc[0];
00261 for (k=0;k<nnsd[j];k++){
00262 mltg[j][k]=buff[k];
00263 }
00264
00265
00266 for (i=1;i<nproc;i++){
00267 MPI_Recv (buff,buffsize,MPI_LONG,MPI_ANY_SOURCE,MPI_ANY_TAG,MPI_COMM_WORLD,&stat);
00268 j=domproc[stat.MPI_TAG];
00269 for (k=0;k<nnsd[j];k++){
00270 mltg[j][k]=buff[k];
00271 }
00272 }
00273 }
00274 else{
00275 MPI_Send (buff,buffsize,MPI_LONG,0,myrank,MPI_COMM_WORLD);
00276 }
00277 MPI_Barrier (MPI_COMM_WORLD);
00278
00279 delete [] buff;
00280
00281
00282 if (myrank==0){
00283 fprintf (out,"\n\n\n mltg (arrays ltg collected on the master)");
00284 for (i=0;i<nproc;i++){
00285 for (j=0;j<nnsd[i];j++){
00286 fprintf (out,"\n mltg %6ld %6ld %ld",i+1,j+1,mltg[i][j]);
00287 }
00288 }
00289 }
00290
00291 }
00292
00293
00294
00295
00296
00297
00298
00299
00300
00301
00302
00303 void partopjk::ltg_conversion (long *ltg,long *domproc,FILE *out,char *proc_name)
00304 {
00305 long i,j,max;
00306
00307 if (myrank==0){
00308
00309
00310
00311
00312 switch (md){
00313 case all_nodes:{
00314 max=0;
00315 for (i=0;i<nproc;i++){
00316 for (j=0;j<nnsd[i];j++){
00317 if (max<mltg[i][j]){
00318 max = mltg[i][j];
00319 }
00320 }
00321 }
00322 break;
00323 }
00324 case bound_nodes:{
00325 max=0;
00326 for (i=0;i<nproc;i++){
00327 for (j=0;j<nnsd[i];j++){
00328 if (max<mltg[i][j]){
00329 max = mltg[i][j];
00330 }
00331 }
00332 }
00333 break;
00334 }
00335 case neg_bound_nodes:{
00336 max=0;
00337 for (i=0;i<nproc;i++){
00338 for (j=0;j<nnsd[i];j++){
00339 if (max<abs(mltg[i][j])){
00340 max = abs(mltg[i][j]);
00341 }
00342 }
00343 }
00344 break;
00345 }
00346 default:{
00347 par_print_err(myrank,proc_name,"wrong type of mesh description is required",__FILE__, __LINE__, __func__);
00348 }
00349 }
00350 max++;
00351
00352
00353
00354
00355 switch (md){
00356 case all_nodes:{
00357
00358 break;
00359 }
00360 case bound_nodes:{
00361 for (i=0;i<nproc;i++){
00362 for (j=0;j<nnsd[i];j++){
00363 if (mltg[i][j]==-1){
00364 mltg[i][j]=max;
00365 max++;
00366 }
00367 }
00368 }
00369 break;
00370 }
00371 case neg_bound_nodes:{
00372 for (i=0;i<nproc;i++){
00373 for (j=0;j<nnsd[i];j++){
00374 if (mltg[i][j]<-1){
00375 mltg[i][j]=0-mltg[i][j];
00376 }
00377 }
00378 }
00379 break;
00380 }
00381 default:{
00382 par_print_err(myrank,proc_name,"wrong type of mesh description is required",__FILE__, __LINE__, __func__);
00383 }
00384 }
00385
00386
00387 tnnp=max;
00388
00389 }
00390
00391
00392 if (myrank==0){
00393 fprintf (out,"\n\n\n tnnp (the total number of nodes in the whole problem) %ld",tnnp);
00394 fprintf (out,"\n\n\n mltg (arrays ltg collected on the master) after conversion");
00395 for (i=0;i<nproc;i++){
00396 for (j=0;j<nnsd[i];j++){
00397 fprintf (out,"\n mltg %6ld %6ld %ld",i+1,j+1,mltg[i][j]);
00398 }
00399 }
00400 }
00401
00402 }
00403
00404
00405
00406
00407
00408
00409
00410
00411 void partopjk::node_multiplicity (FILE *out)
00412 {
00413 long i,j;
00414
00415 if (myrank==0){
00416 nodmultip = new long [tnnp];
00417 for (i=0;i<tnnp;i++){
00418 nodmultip[i]=0;
00419 }
00420
00421 for (i=0;i<nproc;i++){
00422 for (j=0;j<nnsd[i];j++){
00423 nodmultip[mltg[i][j]]++;
00424 }
00425 }
00426 }
00427
00428
00429 if (myrank==0){
00430 fprintf (out,"\n\n\n nodmultip (array of nodal multiplicity for all nodes in the problem)");
00431 for (i=0;i<tnnp;i++){
00432 fprintf (out,"\n nodmultip %6ld %ld",i+1,nodmultip[i]);
00433 }
00434 }
00435
00436 }
00437
00438
00439
00440
00441
00442
00443
00444
00445 void partopjk::node_domain (FILE *out)
00446 {
00447 long i,j,nid;
00448
00449 if (myrank==0){
00450 noddom = new long [tnnp];
00451
00452
00453 for (i=0;i<nproc;i++){
00454
00455 for (j=0;j<nnsd[i];j++){
00456
00457 nid = mltg[i][j];
00458 noddom[nid]=i;
00459 }
00460 }
00461 }
00462
00463
00464 if (myrank==0){
00465 fprintf (out,"\n\n\n noddom (node-domain correspondence)");
00466 for (i=0;i<tnnp;i++){
00467 fprintf (out,"\n noddom %6ld %ld",i+1,noddom[i]);
00468 }
00469 }
00470
00471 }
00472
00473
00474
00475
00476
00477
00478
00479
00480
00481
00482
00483
00484 void partopjk::ndofn_on_master (gtopology *top,long *domproc,FILE *out)
00485 {
00486 long i,j,k,buffsize;
00487 long *buff;
00488 MPI_Status stat;
00489
00490 buffsize = maxnn+1;
00491 buff = new long [buffsize];
00492
00493 maxndofd=0;
00494 for (i=0;i<nn;i++){
00495 buff[i] = top->give_ndofn (i);
00496 maxndofd+=buff[i];
00497 }
00498 buff[nn]=maxndofd;
00499
00500 if (myrank==0){
00501 ndofnm = new long [tnnp];
00502
00503
00504 j=domproc[0];
00505 for (k=0;k<nnsd[j];k++){
00506 ndofnm[mltg[j][k]]=buff[k];
00507 }
00508
00509
00510 for (i=1;i<nproc;i++){
00511 MPI_Recv (buff,buffsize,MPI_LONG,MPI_ANY_SOURCE,MPI_ANY_TAG,MPI_COMM_WORLD,&stat);
00512 j=domproc[stat.MPI_TAG];
00513 for (k=0;k<nnsd[j];k++){
00514 ndofnm[mltg[j][k]]=buff[k];
00515 }
00516 if (maxndofd<buff[nnsd[j]]){
00517 maxndofd=buff[nnsd[j]];
00518 }
00519 }
00520
00521 for (i=1;i<nproc;i++){
00522 MPI_Send (&maxndofd,1,MPI_LONG,i,myrank,MPI_COMM_WORLD);
00523 }
00524 }
00525 else{
00526 MPI_Send (buff,buffsize,MPI_LONG,0,myrank,MPI_COMM_WORLD);
00527 MPI_Recv (&maxndofd,1,MPI_LONG,MPI_ANY_SOURCE,MPI_ANY_TAG,MPI_COMM_WORLD,&stat);
00528 }
00529 MPI_Barrier (MPI_COMM_WORLD);
00530
00531 if (myrank==0){
00532 for (i=1;i<nproc;i++){
00533 MPI_Send (&maxndofd,1,MPI_LONG,i,myrank,MPI_COMM_WORLD);
00534 }
00535 }
00536 else{
00537 MPI_Recv (&maxndofd,1,MPI_LONG,MPI_ANY_SOURCE,MPI_ANY_TAG,MPI_COMM_WORLD,&stat);
00538 }
00539 MPI_Barrier (MPI_COMM_WORLD);
00540
00541 delete [] buff;
00542
00543
00544 fprintf (out,"\n\n maxndofd (maximum number of DOFs on subdomain) %ld",maxndofd);
00545 if (myrank==0){
00546 fprintf (out,"\n\n\n ndofnm (array of nodal ndofn)");
00547 for (i=0;i<tnnp;i++){
00548 fprintf (out,"\n ndofnm %6ld %ld",i+1,ndofnm[i]);
00549 }
00550 }
00551
00552
00553 }
00554
00555
00556
00557
00558
00559
00560
00561
00562
00563
00564 void partopjk::dof_on_master (gtopology *top,long *domproc,FILE *out)
00565 {
00566 long i,j,k,l,m,nid,ndofn,buffsize;
00567 long *buff;
00568 MPI_Status stat;
00569
00570 buffsize=maxndofd;
00571 buff = new long [buffsize];
00572
00573 k=0;
00574 for (i=0;i<nn;i++){
00575 ndofn = top->give_ndofn (i);
00576 for (j=0;j<ndofn;j++){
00577 buff[k]=top->give_dof (i,j);
00578 k++;
00579 }
00580 }
00581
00582 if (myrank==0){
00583 dofm = new long* [tnnp];
00584 for (i=0;i<tnnp;i++){
00585 dofm[i] = new long [ndofnm[i]];
00586 }
00587
00588
00589 j=domproc[0];
00590 l=0;
00591 for (k=0;k<nnsd[j];k++){
00592
00593 nid=mltg[j][k];
00594
00595 ndofn=ndofnm[nid];
00596 for (m=0;m<ndofn;m++){
00597 dofm[nid][m]=buff[l];
00598 l++;
00599 }
00600 }
00601
00602
00603 for (i=1;i<nproc;i++){
00604 MPI_Recv (buff,buffsize,MPI_LONG,MPI_ANY_SOURCE,MPI_ANY_TAG,MPI_COMM_WORLD,&stat);
00605 j=domproc[stat.MPI_TAG];
00606 l=0;
00607 for (k=0;k<nnsd[j];k++){
00608
00609 nid=mltg[j][k];
00610
00611 ndofn=ndofnm[nid];
00612 for (m=0;m<ndofn;m++){
00613 dofm[nid][m]=buff[l];
00614 l++;
00615 }
00616 }
00617 }
00618 }
00619 else{
00620 MPI_Send (buff,buffsize,MPI_LONG,0,myrank,MPI_COMM_WORLD);
00621 }
00622 MPI_Barrier (MPI_COMM_WORLD);
00623
00624 delete [] buff;
00625
00626
00627 if (myrank==0){
00628 fprintf (out,"\n\n\n dofm (array of nodal DOFs on the master)");
00629 for (i=0;i<tnnp;i++){
00630 for (j=0;j<ndofnm[i];j++){
00631 fprintf (out,"\n ndofnm %6ld %6ld %ld",i+1,j+1,dofm[i][j]);
00632 }
00633 }
00634 }
00635
00636
00637 }
00638
00639
00640
00641
00642
00643
00644
00645
00646
00647
00648
00649
00650
00651
00652
00653
00654
00655
00656 void partopjk::coupled_dofs_detection (FILE *out)
00657 {
00658 long i,j,k,max,dofid,stop,nid,nid2,sdid,sdid2;
00659 long *ncdof,**cdof;
00660
00661 if (myrank==0){
00662
00663 max=0;
00664 for (i=0;i<tnnp;i++){
00665 for (j=0;j<ndofnm[i];j++){
00666 if (max<dofm[i][j]){
00667 max=dofm[i][j];
00668 }
00669 }
00670 }
00671
00672 if (max>1){
00673
00674
00675 ncdof = new long [max];
00676 for (i=0;i<max;i++){
00677 ncdof[i]=0;
00678 }
00679 for (i=0;i<tnnp;i++){
00680 for (j=0;j<ndofnm[i];j++){
00681 if (dofm[i][j]>1){
00682 ncdof[dofm[i][j]-2]++;
00683 }
00684 }
00685 }
00686
00687 cdof = new long* [max-1];
00688 for (i=0;i<max-1;i++){
00689 cdof[i] = new long [ncdof[i]];
00690 ncdof[i]=0;
00691 }
00692
00693
00694
00695 for (i=0;i<tnnp;i++){
00696
00697 for (j=0;j<ndofnm[i];j++){
00698 if (dofm[i][j]>1){
00699
00700
00701
00702 dofid=dofm[i][j]-2;
00703
00704 cdof[dofid][ncdof[dofid]]=i;
00705
00706 ncdof[dofid]++;
00707 }
00708 }
00709 }
00710
00711 fprintf (out,"\n\n\n jouda \n\n");
00712 for (i=0;i<max-1;i++){
00713 fprintf (out,"\n ncdof %ld %ld",i,ncdof[i]);
00714 for (j=0;j<ncdof[i];j++){
00715 fprintf (out," cdof %ld %ld %ld",i,j,cdof[i][j]);
00716 }
00717 }
00718 fprintf (out,"\n jouda \n\n\n\n\n");
00719
00720
00721
00722 for (i=0;i<max-1;i++){
00723
00724 stop=0;
00725
00726 for (j=0;j<ncdof[i];j++){
00727
00728 nid=cdof[i][j];
00729 if (nodmultip[nid]>1){
00730
00731
00732
00733
00734
00735 for (k=0;k<ncdof[i];k++){
00736 nid2=cdof[i][k];
00737 if (nodmultip[nid2]==1){
00738 nodmultip[nid2]=-1;
00739 }
00740 }
00741 stop=1;
00742 break;
00743 }
00744 }
00745
00746
00747 if (stop==0){
00748
00749
00750
00751
00752 nid=cdof[i][0];
00753
00754 sdid=noddom[nid];
00755
00756 for (j=0;j<ncdof[i];j++){
00757
00758 nid2=cdof[i][j];
00759
00760 sdid2=noddom[nid2];
00761 if (sdid!=sdid2){
00762
00763
00764
00765 for (k=0;k<ncdof[i];k++){
00766 nid2=cdof[i][k];
00767 if (nodmultip[nid2]==1){
00768 nodmultip[nid2]=-1;
00769 }
00770 }
00771
00772 stop=1;
00773 break;
00774 }
00775
00776 }
00777 }
00778 }
00779
00780
00781 for (i=0;i<max-1;i++){
00782 delete [] cdof[i];
00783 }
00784 delete [] cdof;
00785
00786 delete [] ncdof;
00787
00788 }
00789 }
00790
00791
00792
00793 if (myrank==0){
00794 fprintf (out,"\n\n\n nodmultip (array of nodal multiplicity for all nodes in the problem)");
00795 fprintf (out,"\n the array was modified due to coupled DOFs");
00796 for (i=0;i<tnnp;i++){
00797 fprintf (out,"\n nodmultip %6ld %ld",i+1,nodmultip[i]);
00798 }
00799 }
00800
00801
00802 }
00803
00804
00805
00806
00807
00808
00809
00810
00811 void partopjk::internal_boundary_nodes_detection (FILE *out)
00812 {
00813 long i,j,k,nid,buffsize;
00814 long *buff;
00815 MPI_Status stat;
00816
00817 buffsize = maxnn+2;
00818 buff = new long [buffsize];
00819
00820 if (myrank==0){
00821 fprintf (out,"\n nbnd = %p\n",nbnd);
00822 if (nbnd!=NULL){
00823 delete [] nbnd;
00824 }
00825 nbnd = new long [nproc];
00826 fprintf (out,"\n nbnd = %p\n",nbnd);
00827 if (nind!=NULL){
00828 delete [] nind;
00829 }
00830 nind = new long [nproc];
00831
00832
00833
00834
00835 for (i=1;i<nproc;i++){
00836 nbnd[i]=0;
00837 nind[i]=0;
00838
00839
00840 for (j=0;j<nnsd[i];j++){
00841
00842 nid = mltg[i][j];
00843
00844 if (nodmultip[nid]>1 || nodmultip[nid]<0){
00845 buff[j]=2;
00846 nbnd[i]++;
00847 }
00848 else{
00849 buff[j]=1;
00850 nind[i]++;
00851 }
00852 }
00853 buff[buffsize-2]=nind[i];
00854 buff[buffsize-1]=nbnd[i];
00855
00856 MPI_Send (buff,buffsize,MPI_LONG,i,myrank,MPI_COMM_WORLD);
00857 }
00858
00859
00860 i=0;
00861 nbnd[i]=0;
00862 nind[i]=0;
00863
00864
00865 for (j=0;j<nnsd[i];j++){
00866
00867 nid = mltg[i][j];
00868
00869 if (nodmultip[nid]>1 || nodmultip[nid]<0){
00870 buff[j]=2;
00871 nbnd[i]++;
00872 }
00873 else{
00874 buff[j]=1;
00875 nind[i]++;
00876 }
00877 }
00878 buff[buffsize-2]=nind[i];
00879 buff[buffsize-1]=nbnd[i];
00880 }
00881 else{
00882 MPI_Recv (buff,buffsize,MPI_LONG,MPI_ANY_SOURCE,MPI_ANY_TAG,MPI_COMM_WORLD,&stat);
00883 }
00884 MPI_Barrier (MPI_COMM_WORLD);
00885
00886
00887 nin = buff[buffsize-2];
00888
00889 nbn = buff[buffsize-1];
00890
00891 inid = new long [nin];
00892
00893 bnid = new long [nbn];
00894
00895
00896 j=0; k=0;
00897 for (i=0;i<nn;i++){
00898 if (buff[i]==1){
00899 inid[j]=i;
00900 j++;
00901 }
00902 if (buff[i]==2){
00903 bnid[k]=i;
00904 k++;
00905 }
00906 }
00907
00908 delete [] buff;
00909
00910
00911
00912
00913 if (myrank==0){
00914 fprintf (out,"\n\n\n nbnd (the number of boundary/interface nodes on subdomains)");
00915 fprintf (out,"\n nind (the number of internal nodes on subdomains)",nbn);
00916 for (i=0;i<nproc;i++){
00917 fprintf (out,"\n %ld nbnd %6ld nind %6ld",i+1,nbnd[i],nind[i]);
00918 }
00919 }
00920
00921 fprintf (out,"\n\n\n nin (the number of internal nodes %ld",nin);
00922 fprintf (out,"\n nbn (the number of boundary/interface nodes %ld",nbn);
00923 fprintf (out,"\n\n\n inid (array of internal node numbers)");
00924 for (i=0;i<nin;i++){
00925 fprintf (out,"\n inid %6ld %ld",i+1,inid[i]);
00926 }
00927 fprintf (out,"\n\n\n bnid (array of boundary/interface node numbers)");
00928 for (i=0;i<nbn;i++){
00929 fprintf (out,"\n bnid %6ld %ld",i+1,bnid[i]);
00930 }
00931
00932
00933 }
00934
00935
00936
00937
00938
00939
00940
00941
00942
00943
00944 long partopjk::schur_local_ordering (gtopology *top,FILE *out)
00945 {
00946 long i,j,k,nid,ndofn;
00947 long *aux;
00948
00949
00950
00951
00952 ndof=0;
00953
00954 for (i=0;i<nn;i++){
00955
00956 ndofn = top->give_ndofn (i);
00957
00958 for (j=0;j<ndofn;j++){
00959
00960 k=top->give_dof (i,j);
00961 if (ndof<k){
00962 ndof=k;
00963 }
00964 }
00965 }
00966
00967
00968 ndof--;
00969 if (ndof<0){
00970 ndof=0;
00971 }
00972 aux = new long [ndof];
00973 for (i=0;i<ndof;i++){
00974 aux[i]=-1;
00975 }
00976
00977
00978
00979
00980
00981 ndof=1;
00982
00983 for (i=0;i<nin;i++){
00984
00985 nid=inid[i];
00986
00987 ndofn = top->give_ndofn (nid);
00988
00989 for (j=0;j<ndofn;j++){
00990
00991 k=top->give_dof (nid,j);
00992 if (k<0){
00993
00994 continue;
00995 }
00996 if (k==0){
00997
00998 continue;
00999 }
01000 if (k==1){
01001
01002 top->save_dof (nid,j,ndof);
01003 ndof++;
01004 }
01005 if (k>1){
01006
01007 if (aux[k-2]==-1){
01008 top->save_dof (nid,j,ndof);
01009 aux[k-2]=ndof;
01010 ndof++;
01011 }
01012 else{
01013 top->save_dof (nid,j,aux[k-2]);
01014 }
01015 }
01016 }
01017 }
01018
01019
01020 nidof = ndof-1;
01021
01022
01023
01024
01025
01026 for (i=0;i<nbn;i++){
01027
01028 nid=bnid[i];
01029
01030 ndofn = top->give_ndofn (nid);
01031
01032 for (j=0;j<ndofn;j++){
01033
01034 k=top->give_dof (nid,j);
01035 if (k<0){
01036
01037 continue;
01038 }
01039 if (k==0){
01040
01041 continue;
01042 }
01043 if (k==1){
01044
01045 top->save_dof (nid,j,ndof);
01046 ndof++;
01047 }
01048 if (k>1){
01049
01050 if (aux[k-2]==-1){
01051 top->save_dof (nid,j,ndof);
01052 aux[k-2]=ndof;
01053 ndof++;
01054 }
01055 else{
01056 top->save_dof (nid,j,aux[k-2]);
01057 }
01058 }
01059 }
01060 }
01061
01062 delete [] aux;
01063
01064
01065 ndof--;
01066
01067
01068 nbdof = ndof-nidof;
01069
01070
01071
01072 top->cnstate=1;
01073
01074
01075 fprintf (out,"\n\n\n nidof (the number of internal DOFs %ld",nidof);
01076 fprintf (out,"\n nbdof (the number of boundary/interface DOFs %ld",nbdof);
01077 fprintf (out,"\n ndof (the number of all DOFs %ld",ndof);
01078 fprintf (out,"\n\n\n code number");
01079
01080 for (i=0;i<nn;i++){
01081 fprintf (out,"\n node %6ld ",i+1);
01082
01083 ndofn = top->give_ndofn (i);
01084
01085 for (j=0;j<ndofn;j++){
01086 fprintf (out," %6ld",top->give_dof (i,j));
01087 }
01088 }
01089 fprintf (out,"\n\n");
01090
01091
01092 return ndof;
01093 }
01094
01095
01096
01097
01098
01099
01100
01101
01102
01103
01104
01105 void partopjk::schur_coarse_ordering (FILE *out)
01106 {
01107 long i,j,k;
01108 long *aux;
01109
01110 if (myrank==0){
01111
01112 ndofc=0;
01113
01114
01115 for (i=0;i<tnnp;i++){
01116 if (nodmultip[i]>1 || nodmultip[i]<0){
01117
01118
01119
01120
01121 for (j=0;j<ndofnm[i];j++){
01122 if (ndofc<dofm[i][j]){
01123 ndofc=dofm[i][j];
01124 }
01125 }
01126 }
01127 }
01128
01129
01130 ndofc--;
01131 if (ndofc<0){
01132 ndofc=0;
01133 }
01134 aux = new long [ndofc];
01135 for (i=0;i<ndofc;i++){
01136 aux[i]=-1;
01137 }
01138
01139 ndofc=1;
01140
01141 for (i=0;i<tnnp;i++){
01142 if (nodmultip[i]>1 || nodmultip[i]<0){
01143
01144
01145
01146
01147 for (j=0;j<ndofnm[i];j++){
01148
01149 k=dofm[i][j];
01150 if (k<0){
01151
01152 continue;
01153 }
01154 if (k==0){
01155
01156 continue;
01157 }
01158 if (k==1){
01159
01160 dofm[i][j]=ndofc;
01161 ndofc++;
01162 }
01163 if (k>1){
01164
01165 if (aux[k-2]==-1){
01166 dofm[i][j]=ndofc;
01167 aux[k-2]=ndofc;
01168 ndofc++;
01169 }
01170 else{
01171 dofm[i][j]=aux[k-2];
01172 }
01173 }
01174 }
01175 }
01176 }
01177 ndofc--;
01178
01179 delete [] aux;
01180
01181 }
01182
01183
01184
01185
01186
01187 for (i=0;i<tnnp;i++){
01188 if (nodmultip[i]==1){
01189 for (j=0;j<ndofnm[i];j++){
01190 dofm[i][j]=0;
01191 }
01192 }
01193 }
01194
01195
01196 if (myrank==0){
01197 fprintf (out,"\n\n\n ndofc %ld",ndofc);
01198 fprintf (out,"\n\n\n dofm after code numbers generation ");
01199 fprintf (out,"\n dofm (array of nodal DOFs on the master)");
01200 for (i=0;i<tnnp;i++){
01201 for (j=0;j<ndofnm[i];j++){
01202 fprintf (out,"\n dofm %6ld %6ld %ld",i+1,j+1,dofm[i][j]);
01203 }
01204 }
01205 }
01206
01207
01208 }
01209
01210
01211
01212
01213
01214
01215
01216
01217
01218
01219
01220 void partopjk::coarse_code_numbers (FILE *out)
01221 {
01222 long i,j,k,nid,dofid,buffsize;
01223 long *aux,*buff;
01224 MPI_Status stat;
01225
01226 if (myrank==0){
01227 if (nbdofd!=NULL){
01228 delete [] nbdofd;
01229 }
01230 nbdofd = new long [nproc];
01231
01232 aux = new long [ndofc];
01233
01234
01235 for (i=0;i<nproc;i++){
01236 for (j=0;j<ndofc;j++){
01237 aux[j]=0;
01238 }
01239
01240 nbdofd[i]=0;
01241
01242
01243 for (j=0;j<nnsd[i];j++){
01244
01245 nid = mltg[i][j];
01246 if (nodmultip[nid]>1 || nodmultip[nid]<0){
01247
01248
01249 for (k=0;k<ndofnm[nid];k++){
01250
01251 dofid=dofm[nid][k];
01252 if (dofid>0){
01253 aux[dofid-1]++;
01254 if (aux[dofid-1]==1){
01255 nbdofd[i]++;
01256 }
01257 }
01258 }
01259 }
01260 }
01261 }
01262
01263 if (bdofd!=NULL){
01264 for (i=0;i<nproc;i++){
01265 delete [] bdofd[i];
01266 }
01267 delete [] bdofd;
01268 }
01269 bdofd = new long* [nproc];
01270 for (i=0;i<nproc;i++){
01271 fprintf (stdout,"\n nbdofd %ld \n",nbdofd[i]);
01272 bdofd[i] = new long [nbdofd[i]];
01273 }
01274
01275
01276
01277 for (i=0;i<nproc;i++){
01278 for (j=0;j<ndofc;j++){
01279 aux[j]=0;
01280 }
01281
01282 nbdofd[i]=0;
01283
01284
01285 for (j=0;j<nnsd[i];j++){
01286
01287 nid = mltg[i][j];
01288 if (nodmultip[nid]>1 || nodmultip[nid]<0){
01289
01290
01291 for (k=0;k<ndofnm[nid];k++){
01292
01293 dofid=dofm[nid][k];
01294 if (dofid>0){
01295 aux[dofid-1]++;
01296 if (aux[dofid-1]==1){
01297 bdofd[i][nbdofd[i]]=dofm[nid][k];
01298 nbdofd[i]++;
01299 }
01300 }
01301 }
01302 }
01303 }
01304 }
01305
01306 delete [] aux;
01307 }
01308
01309
01310 if (myrank==0){
01311 maxnbdofd = 0;
01312 for (i=0;i<nproc;i++){
01313 fprintf (stdout,"\n nbdofd %ld \n",nbdofd[i]);
01314 if (maxnbdofd<nbdofd[i]){
01315 maxnbdofd=nbdofd[i];
01316 }
01317 }
01318 for (i=1;i<nproc;i++){
01319 MPI_Send (&maxnbdofd,1,MPI_LONG,i,myrank,MPI_COMM_WORLD);
01320 }
01321 }
01322 else{
01323 MPI_Recv (&maxnbdofd,1,MPI_LONG,MPI_ANY_SOURCE,MPI_ANY_TAG,MPI_COMM_WORLD,&stat);
01324 }
01325 MPI_Barrier (MPI_COMM_WORLD);
01326
01327
01328
01329 fprintf (out,"\n\n\n maxnbdofd (numbers of boundary/interface DOFs on subdomains) %ld",maxnbdofd);
01330 if (myrank==0){
01331 fprintf (out,"\n\n\n nbdofd (numbers of boundary/interface DOFs on subdomains)");
01332 for (i=0;i<nproc;i++){
01333 fprintf (out,"\n nbdofd %6ld %ld",i+1,nbdofd[i]);
01334 }
01335 fprintf (out,"\n\n\n nbdofd (numbers of boundary/interface DOFs on subdomains)");
01336 for (i=0;i<nproc;i++){
01337 for (j=0;j<nbdofd[i];j++){
01338 fprintf (out,"\n bdofd %6ld %6ld %ld",i+1,j+1,bdofd[i][j]);
01339 }
01340 }
01341 }
01342
01343
01344
01345 buffsize=maxnbdofd;
01346 buff = new long [buffsize];
01347
01348 if (myrank==0){
01349
01350 fprintf (out,"\n nbnd = %p\n",nbnd);
01351 for (i=1;i<nproc;i++){
01352 fprintf (out," smycka i %ld",i);
01353 for (j=0;j<nbdofd[i];j++){
01354 buff[j]=bdofd[i][j];
01355 }
01356 MPI_Send (buff,buffsize,MPI_LONG,i,myrank,MPI_COMM_WORLD);
01357 }
01358
01359 i=0;
01360 for (j=0;j<nbdofd[i];j++){
01361 buff[j]=bdofd[i][j];
01362 }
01363
01364 }
01365 else{
01366 MPI_Recv (buff,buffsize,MPI_LONG,MPI_ANY_SOURCE,MPI_ANY_TAG,MPI_COMM_WORLD,&stat);
01367 }
01368 MPI_Barrier (MPI_COMM_WORLD);
01369
01370
01371 bdof = new long [nbdof];
01372 for (i=0;i<nbdof;i++){
01373 bdof[i]=buff[i];
01374 }
01375
01376 delete [] buff;
01377
01378
01379
01380 fprintf (out,"\n\n\n bdof (boundary/interface DOFs on subdomain");
01381 for (j=0;j<nbdof;j++){
01382 fprintf (out,"\n bdof %6ld %ld",j+1,bdof[j]);
01383 }
01384
01385
01386
01387 }
01388
01389
01390
01391 long partopjk::vse (long *ltg,gtopology *top,long *domproc,FILE *out,char *proc_name)
01392 {
01393 long ndof;
01394
01395 fprintf (out,"\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n ZACATEK NOVEHO PARTOPU \n\n\n\n\n\n\n\n\n");
01396
01397
01398 numbers_of_nodes_on_subdomains (top,domproc,out);
01399
01400
01401 ltg_on_master (ltg,domproc,out);
01402
01403
01404 ltg_conversion (ltg,domproc,out,proc_name);
01405
01406
01407 node_multiplicity (out);
01408
01409
01410 node_domain (out);
01411
01412
01413 ndofn_on_master (top,domproc,out);
01414
01415
01416 dof_on_master (top,domproc,out);
01417
01418
01419 coupled_dofs_detection (out);
01420
01421
01422 internal_boundary_nodes_detection (out);
01423
01424
01425 ndof = schur_local_ordering (top,out);
01426
01427
01428 schur_coarse_ordering (out);
01429
01430
01431 coarse_code_numbers (out);
01432
01433 fprintf (out,"\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n KONEC NOVEHO PARTOPU \n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n");
01434
01435 return ndof;
01436
01437 }
01438
01439
01440
01441
01442
01443
01444
01445
01446
01447
01448
01449
01450
01451 void partopjk::full_node_domain (FILE *out)
01452 {
01453 long i,j,k,l,min;
01454
01455 if (myrank==0){
01456 if (fullnoddom!=NULL){
01457 for (i=0;i<tnnp;i++){
01458 delete [] fullnoddom[i];
01459 }
01460 delete [] fullnoddom;
01461 }
01462 fullnoddom = new long* [tnnp];
01463 for (i=0;i<tnnp;i++){
01464 fullnoddom[i] = new long [nodmultip[i]];
01465
01466
01467 nodmultip[i]=0;
01468 }
01469
01470
01471 for (i=0;i<nproc;i++){
01472
01473 for (j=0;j<nnsd[i];j++){
01474
01475 k=mltg[i][j];
01476
01477 l=nodmultip[k];
01478 fullnoddom[k][l]=i;
01479 nodmultip[k]++;
01480 }
01481 }
01482
01483
01484
01485 for (i=0;i<tnnp;i++){
01486
01487 for (j=0;j<nodmultip[i];j++){
01488 min=nproc; l=-1;
01489 for (k=j;k<nodmultip[i];k++){
01490 if (min>fullnoddom[i][k]){
01491 min=fullnoddom[i][k];
01492 l=k;
01493 }
01494 }
01495 min=fullnoddom[i][j];
01496 fullnoddom[i][j]=fullnoddom[i][l];
01497 fullnoddom[i][l]=min;
01498 }
01499 }
01500
01501 }
01502
01503
01504
01505 if (myrank==0){
01506 fprintf (out,"\n\n\n fullnoddom (array of full node-domain correspondence)");
01507 for (i=0;i<tnnp;i++){
01508 for (j=0;j<nodmultip[i];j++){
01509 fprintf (out,"\n fullnoddom %6ld %6ld %ld",i+1,j+1,fullnoddom[i][j]);
01510 }
01511 }
01512 }
01513
01514
01515 }
01516
01517
01518
01519
01520
01521
01522
01523
01524
01525
01526
01527 void partopjk::lagrange_multipliers (psolver ps,FILE *out)
01528 {
01529 long i,j,k,l,m,ndofn;
01530
01531 if (myrank==0){
01532 lagrmultip = new long* [tnnp];
01533 for (i=0;i<tnnp;i++){
01534 if (ps.f1->fetiimpl == nonredundant){
01535 lagrmultip[i] = new long [(nodmultip[i]-1)*ndofnm[i]];
01536 for (j=0;j<(nodmultip[i]-1)*ndofnm[i];j++){
01537 lagrmultip[i][j]=0;
01538 }
01539 }
01540 if (ps.f1->fetiimpl == redundant){
01541 lagrmultip[i] = new long [nodmultip[i]*(nodmultip[i]-1)*ndofnm[i]/2];
01542 for (j=0;j<nodmultip[i]*(nodmultip[i]-1)*ndofnm[i]/2;j++){
01543 lagrmultip[i][j]=0;
01544 }
01545 }
01546 }
01547
01548
01549 ndofc=1;
01550
01551
01552 for (i=0;i<tnnp;i++){
01553 if (nodmultip[i]>1){
01554
01555
01556 if (ps.f1->fetiimpl == nonredundant){
01557
01558 ndofn=ndofnm[i];
01559 for (j=0;j<ndofn;j++){
01560 if (dofm[i][j]>0){
01561
01562 for (k=0;k<nodmultip[i]-1;k++){
01563 lagrmultip[i][k*ndofn+j]=ndofc;
01564 ndofc++;
01565 }
01566 }
01567 }
01568 }
01569
01570 if (ps.f1->fetiimpl == redundant){
01571
01572 ndofn=ndofnm[i];
01573 for (j=0;j<ndofn;j++){
01574 if (dofm[i][j]>0){
01575 m=0;
01576
01577
01578 for (k=0;k<nodmultip[i]-1;k++){
01579 for (l=k;l<nodmultip[i]-1;l++){
01580 lagrmultip[i][m*ndofn+j]=ndofc;
01581 ndofc++;
01582 m++;
01583 }
01584 }
01585 }
01586 }
01587 }
01588
01589
01590 }
01591 }
01592
01593 }
01594
01595
01596
01597 if (myrank==0){
01598 fprintf (out,"\n\n\n lagrmultip ()");
01599 for (i=0;i<tnnp;i++){
01600 fprintf (out,"\n node %8ld (nodmultip %2ld)",i+1,nodmultip[i]);
01601
01602 if (ps.f1->fetiimpl == nonredundant){
01603 for (j=0;j<(nodmultip[i]-1)*ndofnm[i];j++){
01604 fprintf (out," %ld",lagrmultip[i][j]);
01605 }
01606 }
01607
01608 if (ps.f1->fetiimpl == redundant){
01609 for (j=0;j<nodmultip[i]*(nodmultip[i]-1)*ndofnm[i]/2;j++){
01610 fprintf (out," %ld",lagrmultip[i][j]);
01611 }
01612 }
01613
01614 }
01615 fprintf (out,"\n\n");
01616 }
01617
01618
01619 }
01620
01621
01622
01623
01624
01625
01626
01627
01628
01629
01630
01631
01632
01633 void partopjk::coarse_code_numbers_feti (psolver ps,FILE *out)
01634 {
01635 long i,j,k,l,nid,aux,buffsize,ndofn;
01636 long *buff;
01637 MPI_Status stat;
01638
01639 if (myrank==0){
01640 maxnbdofd=0;
01641
01642 for (i=0;i<nproc;i++){
01643 aux=0;
01644
01645 for (j=0;j<nnsd[i];j++){
01646
01647 nid=mltg[i][j];
01648 if (nodmultip[nid]>1){
01649
01650 if (ps.f1->fetiimpl == nonredundant){
01651 if (fullnoddom[nid][0]==i){
01652 for (k=0;k<ndofnm[nid];k++){
01653 if (dofm[nid][k]>0){
01654 aux+=nodmultip[i]-1;
01655 }
01656 }
01657 }
01658 else{
01659 for (k=0;k<ndofnm[nid];k++){
01660 if (dofm[nid][k]>0){
01661 aux++;
01662 }
01663 }
01664 }
01665 }
01666
01667 if (ps.f1->fetiimpl == redundant){
01668 for (k=0;k<ndofnm[nid];k++){
01669 if (dofm[nid][k]>0){
01670 aux+=nodmultip[i]-1;
01671 }
01672 }
01673 }
01674
01675 }
01676 }
01677
01678 if (maxnbdofd<aux){
01679 maxnbdofd=aux;
01680 }
01681 }
01682 }
01683
01684
01685 if (myrank==0){
01686
01687 for (i=1;i<nproc;i++){
01688 MPI_Send (&maxnbdofd,1,MPI_LONG,i,myrank,MPI_COMM_WORLD);
01689 }
01690 }
01691 else{
01692 MPI_Recv (&maxnbdofd,1,MPI_LONG,MPI_ANY_SOURCE,MPI_ANY_TAG,MPI_COMM_WORLD,&stat);
01693 }
01694 MPI_Barrier (MPI_COMM_WORLD);
01695
01696
01697 buffsize=maxnbdofd;
01698 buff=new long [buffsize];
01699
01700
01701 if (myrank==0){
01702
01703 for (i=1;i<nproc;i++){
01704 aux=0;
01705
01706 for (j=0;j<nnsd[i];j++){
01707
01708 nid=mltg[i][j];
01709 if (nodmultip[nid]>1){
01710
01711 if (ps.f1->fetiimpl == nonredundant){
01712 if (fullnoddom[nid][0]==i){
01713 ndofn=ndofnm[nid];
01714 for (k=0;k<ndofn;k++){
01715 if (dofm[nid][k]>0){
01716 for (l=0;l<nodmultip[nid];l++){
01717 buff[aux]=lagrmultip[nid][l*ndofn+k];
01718 aux++;
01719 }
01720 }
01721 }
01722 }
01723 else{
01724 for (k=0;k<ndofnm[nid];k++){
01725 if (dofm[nid][k]>0){
01726 for (l=0;l<nodmultip[nid];l++){
01727 if (fullnoddom[nid][l]==i){
01728 break;
01729 }
01730 }
01731 buff[aux]=lagrmultip[nid][(l-1)*ndofn+k];
01732 aux++;
01733 }
01734 }
01735 }
01736 }
01737
01738 if (ps.f1->fetiimpl == redundant){
01739 for (k=0;k<ndofnm[nid];k++){
01740 if (dofm[nid][k]>0){
01741 aux+=nodmultip[i]-1;
01742 }
01743 }
01744 }
01745
01746 }
01747 }
01748
01749 if (maxnbdofd<aux){
01750 maxnbdofd=aux;
01751 }
01752 }
01753 }
01754
01755
01756
01757
01758 }