00001 #include "lplate.h"
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012 lplate::lplate (int np,int mr,int nd,gtopology *top)
00013 {
00014 long i,j,k;
00015
00016 zero=1.0e-12;
00017
00018 nproc=np; myrank=mr; ndom=nd;
00019
00020 nl = nproc;
00021 nn = top->nn;
00022 ndofn = top->gnodes[0].ndofn;
00023 nnmult = 4;
00024
00025 if (myrank==0){
00026 lgnodes = new lgnode [nn];
00027 for (i=0;i<nn;i++){
00028 lgnodes[i].cn = new long* [nl+1];
00029 for (j=0;j<=nl;j++){
00030 lgnodes[i].cn[j] = new long [nnmult];
00031 for (k=0;k<nnmult;k++){
00032 lgnodes[i].cn[j][k]=0;
00033 }
00034 }
00035 }
00036
00037 cn = new long** [nl];
00038 for (i=0;i<nl;i++){
00039 cn[i] = new long* [nn];
00040 for (j=0;j<nn;j++){
00041 cn[i][j] = new long [ndofn];
00042 for (k=0;k<ndofn;k++){
00043 cn[i][j][k]=0;
00044 }
00045 }
00046 }
00047
00048 thick = new double* [nl];
00049 for (i=0;i<nl;i++){
00050 thick[i] = new double [nn];
00051 }
00052
00053 constrmat = new double* [(nl-1)*nn*nnmult];
00054 for (i=0;i<(nl-1)*nn*nnmult;i++){
00055 constrmat[i] = new double [nl*2];
00056 for (j=0;j<nl*2;j++){
00057 constrmat[i][j]=0.0;
00058 }
00059 }
00060
00061 }
00062
00063 }
00064
00065 lplate::~lplate ()
00066 {
00067 long i,j;
00068
00069 if (myrank==0){
00070
00071 for (i=0;i<nl;i++){
00072 delete [] thick[i];
00073 }
00074 delete [] thick;
00075
00076 for (i=0;i<nl;i++){
00077 for (j=0;j<nn;j++){
00078 delete [] cn[i][j];
00079 }
00080 delete [] cn[i];
00081 }
00082 delete [] cn;
00083
00084 for (i=0;i<(nl-1)*nn*nnmult;i++){
00085 delete [] constrmat[i];
00086 }
00087 delete [] constrmat;
00088 }
00089
00090 }
00091
00092
00093
00094
00095
00096
00097
00098
00099
00100
00101
00102
00103
00104
00105
00106
00107 void lplate::globcnnum_lpp (gtopology *top,long *domproc,FILE *out)
00108 {
00109 long i,j,k,buffsize;
00110 long *buff;
00111 MPI_Status stat;
00112
00113
00114 buffsize = nn*ndofn;
00115 buff = new long [buffsize];
00116
00117 k=0;
00118 for (i=0;i<nn;i++){
00119 for (j=0;j<ndofn;j++){
00120 buff[k]=top->gnodes[i].cn[j];
00121 k++;
00122 }
00123 }
00124
00125 if (myrank==0){
00126 assemble_cn (buff,ndom);
00127
00128 for (i=1;i<nproc;i++){
00129 MPI_Recv (buff,buffsize,MPI_LONG,MPI_ANY_SOURCE,MPI_ANY_TAG,MPI_COMM_WORLD,&stat);
00130 assemble_cn (buff,domproc[stat.MPI_TAG]);
00131 }
00132
00133
00134 maxndofdom ();
00135
00136 for (i=1;i<nproc;i++){
00137 MPI_Send (&maxndof,1,MPI_LONG,i,myrank,MPI_COMM_WORLD);
00138 }
00139
00140
00141
00142
00143
00144
00145
00146
00147
00148
00149
00150
00151
00152
00153
00154
00155
00156 codenum_mult (nmult);
00157
00158
00159
00160
00161
00162
00163
00164
00165
00166
00167
00168
00169
00170
00171
00172
00173
00174
00175
00176
00177
00178
00179
00180
00181
00182
00183
00184
00185
00186
00187
00188
00189
00190
00191
00192
00193
00194
00195 }
00196 else{
00197 MPI_Send (buff,buffsize,MPI_LONG,0,myrank,MPI_COMM_WORLD);
00198 MPI_Recv (&maxndof,1,MPI_LONG,MPI_ANY_SOURCE,MPI_ANY_TAG,MPI_COMM_WORLD,&stat);
00199 }
00200 MPI_Barrier (MPI_COMM_WORLD);
00201 delete [] buff;
00202
00203 }
00204
00205 void lplate::globcnn_auxiliary ()
00206 {
00207 long i,j,k;
00208
00209 gcn = new long** [nl];
00210 for (i=0;i<nl;i++){
00211 gcn[i] = new long* [nn];
00212 for (j=0;j<nn;j++){
00213 gcn[i][j] = new long [ndofn];
00214 for (k=0;k<ndofn;k++){
00215 gcn[i][j][k]=0;
00216 }
00217 }
00218 }
00219
00220 adr = new long [nproc+1];
00221 adr[0]=0;
00222
00223 ngdof=1;
00224 for (i=0;i<nl;i++){
00225 for (j=0;j<nn;j++){
00226 for (k=0;k<ndofn;k++){
00227 if (cn[i][j][k]>0){
00228 gcn[i][j][k]=ngdof;
00229 ngdof++;
00230 }
00231 }
00232 }
00233 adr[i+1]=ngdof-1;
00234 }
00235 ngdof--;
00236
00237 }
00238
00239
00240
00241
00242
00243
00244
00245
00246
00247
00248
00249
00250
00251 void lplate::constraintmat (double *th,long *domproc,FILE *out)
00252 {
00253 long i;
00254 double *rbuff;
00255 MPI_Status stat;
00256
00257
00258 rbuff = new double [nn];
00259 for (i=0;i<nn;i++){
00260 rbuff[i]=th[i];
00261 }
00262
00263 if (myrank==0){
00264 assemble_thicknesses (rbuff,ndom);
00265
00266 for (i=1;i<nproc;i++){
00267 MPI_Recv (rbuff,nn,MPI_DOUBLE,MPI_ANY_SOURCE,MPI_ANY_TAG,MPI_COMM_WORLD,&stat);
00268 assemble_thicknesses (rbuff,domproc[stat.MPI_TAG]);
00269 }
00270 }
00271 else{
00272 MPI_Send (rbuff,nn,MPI_DOUBLE,0,myrank,MPI_COMM_WORLD);
00273 }
00274 MPI_Barrier (MPI_COMM_WORLD);
00275 delete [] rbuff;
00276
00277
00278
00279 if (myrank==0){
00280 assemble_constr ();
00281
00282
00283
00284 }
00285 }
00286
00287
00288
00289
00290
00291
00292
00293 void lplate::maxndofdom ()
00294 {
00295 long i,j,k;
00296
00297 maxndof=0;
00298 for (i=0;i<nproc;i++){
00299 for (j=0;j<nn;j++){
00300 for (k=0;k<ndofn;k++){
00301 if (cn[i][j][k]>maxndof) maxndof=cn[i][j][k];
00302 }
00303 }
00304 }
00305 }
00306
00307
00308
00309
00310
00311
00312
00313
00314
00315
00316 void lplate::assemble_cn (long *buff,long ndom)
00317 {
00318 long i,j,k;
00319
00320 k=0;
00321 for (i=0;i<nn;i++){
00322 for (j=0;j<ndofn;j++){
00323 cn[ndom][i][j]=buff[k];
00324 k++;
00325 }
00326 }
00327 }
00328
00329
00330
00331
00332
00333
00334
00335
00336
00337
00338 void lplate::assemble_thicknesses (double *buff,long ndom)
00339 {
00340 long i;
00341
00342 for (i=0;i<nn;i++){
00343 thick[ndom][i]=buff[i];
00344 }
00345 }
00346
00347
00348
00349
00350
00351
00352
00353
00354
00355 void lplate::codenum_mult (long &nm)
00356 {
00357 long i,j;
00358
00359 nm=1;
00360 for (i=1;i<nl;i++){
00361 for (j=0;j<nn;j++){
00362 if (cn[i-1][j][0]>0 || cn[i-1][j][4]>0 || cn[i][j][0]>0 || cn[i][j][4]>0){
00363 lgnodes[j].cn[i][0]=nm; nm++;
00364 }
00365 if (cn[i-1][j][1]>0 || cn[i-1][j][3]>0 || cn[i][j][1]>0 || cn[i][j][3]>0){
00366 lgnodes[j].cn[i][1]=nm; nm++;
00367 }
00368 if (cn[i-1][j][2]>0 || cn[i][j][2]>0){
00369 lgnodes[j].cn[i][2]=nm; nm++;
00370 }
00371 if (cn[i-1][j][5]>0 || cn[i][j][5]>0){
00372 lgnodes[j].cn[i][3]=nm; nm++;
00373 }
00374 }
00375 }
00376 nm--;
00377 }
00378
00379
00380
00381
00382
00383
00384
00385
00386
00387
00388
00389
00390
00391
00392
00393 void lplate::assemble_constr ()
00394 {
00395 long i,j,ii;
00396 double lthickness,uthickness;
00397
00398 ii=0;
00399 for (i=1;i<nl;i++){
00400 for (j=0;j<nn;j++){
00401 lthickness = thick[i-1][j];
00402 uthickness = thick[i][j];
00403
00404 if (cn[i-1][j][0]>0) constrmat[ii][(i-1)*2+0]=-1.0;
00405 if (cn[i-1][j][4]>0) constrmat[ii][(i-1)*2+1]=lthickness/(-2.0);
00406 if (cn[i][j][0]>0) constrmat[ii][i*2+0]=1.0;
00407 if (cn[i][j][4]>0) constrmat[ii][i*2+1]=uthickness/(-2.0);
00408 ii++;
00409
00410 if (cn[i-1][j][1]>0) constrmat[ii][(i-1)*2+0]=-1.0;
00411 if (cn[i-1][j][3]>0) constrmat[ii][(i-1)*2+1]=lthickness/2.0;
00412 if (cn[i][j][1]>0) constrmat[ii][i*2+0]=1.0;
00413 if (cn[i][j][3]>0) constrmat[ii][i*2+1]=uthickness/2.0;
00414 ii++;
00415
00416 if (cn[i-1][j][2]>0) constrmat[ii][(i-1)*2+0]=-1.0;
00417 if (cn[i][j][2]>0) constrmat[ii][i*2+0]=1.0;
00418 ii++;
00419
00420 if (cn[i-1][j][5]>0) constrmat[ii][(i-1)*2+0]=-1.0;
00421 if (cn[i][j][5]>0) constrmat[ii][i*2+0]=1.0;
00422 ii++;
00423 }
00424 }
00425
00426 }
00427
00428
00429
00430
00431
00432
00433 void lplate::assemble_constr_dm ()
00434 {
00435 long i,j,k,l,*cnd,*cnm;
00436 double tl,tu;
00437 matrix lcm;
00438
00439 cnd = new long [2*ndofn];
00440 cnm = new long [nnmult];
00441 allocm (nnmult,2*ndofn,lcm);
00442
00443 fprintf (stdout,"\n pocet vsech multiplikatoru %ld",nmult);
00444 fprintf (stdout,"\n pocet vsech posunu %ld",ndof);
00445 fprintf (stdout,"\n pocet vsech posunu %ld",ndof*nproc);
00446
00447
00448 allocm (nmult,ngdof,dcm);
00449
00450 for (i=1;i<nl;i++){
00451 for (j=0;j<nn;j++){
00452
00453 tl = thick[i-1][j];
00454 tu = thick[i][j];
00455
00456 for (k=0;k<ndofn;k++){
00457 cnd[k]=gcn[i-1][j][k];
00458 }
00459 for (k=0;k<ndofn;k++){
00460 cnd[k+ndofn]=gcn[i][j][k];
00461 }
00462 for (k=0;k<nnmult;k++){
00463 cnm[k]=lgnodes[j].cn[i][k];
00464 }
00465
00466
00467
00468
00469
00470
00471
00472
00473
00474
00475
00476
00477
00478 fillm (0.0,lcm);
00479
00480 lcm[0][0] = -1.0;
00481 lcm[1][1] = -1.0;
00482 lcm[2][2] = -1.0;
00483 lcm[1][3] = tl/2.0;
00484 lcm[0][4] = tl/(-2.0);
00485 lcm[3][5] = -1.0;
00486
00487 lcm[0][6] = 1.0;
00488 lcm[1][7] = 1.0;
00489 lcm[2][8] = 1.0;
00490 lcm[1][9] = tu/2.0;
00491 lcm[0][10] = tu/(-2.0);
00492 lcm[3][11] = 1.0;
00493
00494
00495 for (k=0;k<nnmult;k++){
00496 if (cnm[k]==0) continue;
00497 for (l=0;l<2*ndofn;l++){
00498 if (cnd[l]==0) continue;
00499 else dcm[cnm[k]-1][cnd[l]-1]=lcm[k][l];
00500 }
00501 }
00502
00503 }
00504 }
00505
00506 delete [] cnd;
00507 delete [] cnm;
00508 destrm (lcm);
00509
00510 }
00511
00512
00513
00514
00515
00516
00517
00518
00519
00520
00521
00522
00523 void lplate::orthonormalization ()
00524 {
00525 long i,j,k,ii,jj;
00526 double alpha,nor;
00527
00528
00529 for (i=0;i<nn*nnmult;i++){
00530 nor = constrmat[i][0]*constrmat[i][0] + constrmat[i][1]*constrmat[i][1];
00531 nor += constrmat[i][2]*constrmat[i][2] + constrmat[i][3]*constrmat[i][3];
00532 if (nor>zero){
00533 nor=sqrt(nor);
00534 constrmat[i][0]/=nor;
00535 constrmat[i][1]/=nor;
00536 constrmat[i][2]/=nor;
00537 constrmat[i][3]/=nor;
00538 }
00539 }
00540
00541 ii=nn*nnmult;
00542 for (i=1;i<nl-1;i++){
00543 jj=ii-nn*nnmult;
00544 for (j=0;j<nn*nnmult;j++){
00545 alpha=constrmat[ii][i*2+0]*constrmat[jj][i*2+0] + constrmat[ii][i*2+1]*constrmat[jj][i*2+1];
00546 nor=0.0;
00547 for (k=0;k<2*(i+2);k++){
00548 constrmat[ii][k]-=alpha*constrmat[jj][k];
00549 nor+=constrmat[ii][k]*constrmat[ii][k];
00550 }
00551 if (nor>zero){
00552 nor=sqrt(nor);
00553 for (k=0;k<2*(i+2);k++){
00554 constrmat[ii][k]/=nor;
00555 }
00556 }
00557 ii++; jj++;
00558 }
00559 }
00560 }
00561
00562
00563
00564
00565
00566
00567 void lplate::orthonormalization_dm ()
00568 {
00569 long i,j,k;
00570 double nor,alpha;
00571
00572 nor=0.0;
00573 for (i=0;i<ngdof;i++){
00574 nor+=dcm[0][i]*dcm[0][i];
00575 }
00576 nor=sqrt(nor);
00577 if (nor>zero){
00578 for (i=0;i<ngdof;i++){
00579 dcm[0][i]/=nor;
00580 }
00581 }
00582
00583 for (i=0;i<nmult;i++){
00584 for (j=i+1;j<nmult;j++){
00585 alpha=0.0;
00586 for (k=0;k<ngdof;k++){
00587 alpha+=dcm[j][k]*dcm[i][k];
00588 }
00589 nor=0.0;
00590 for (k=0;k<ngdof;k++){
00591 dcm[j][k]-=alpha*dcm[i][k];
00592 nor+=dcm[j][k]*dcm[j][k];
00593 }
00594 nor=sqrt(nor);
00595 if (nor>zero){
00596 for (k=0;k<ngdof;k++){
00597 dcm[j][k]/=nor;
00598 }
00599 }
00600 }
00601 }
00602
00603 }
00604
00605
00606
00607
00608
00609
00610
00611
00612
00613
00614
00615
00616
00617 void lplate::locglob (double *gv,double *lv,long ns)
00618 {
00619 long i,j,li,ii,cnl,cn1,cn2;
00620 double cm1,cm2,lv1,lv2;
00621
00622 li=ns-1;
00623 if (li<0) li=0;
00624 ii=li*nn*nnmult;
00625 for (i=li;i<nl-1;i++){
00626 for (j=0;j<nn;j++){
00627
00628 cnl=lgnodes[j].cn[i+1][0];
00629 if (cnl==0){ ii++; }
00630 else{
00631 cm1=constrmat[ii][ns*2+0];
00632 cm2=constrmat[ii][ns*2+1];
00633 cn1=cn[ns][j][0];
00634 if (cn1>0) lv1=lv[cn1-1];
00635 else lv1=0.0;
00636 cn2=cn[ns][j][4];
00637 if (cn2>0) lv2=lv[cn2-1];
00638 else lv2=0.0;
00639 gv[cnl-1]+=cm1*lv1+cm2*lv2;
00640 ii++;
00641 }
00642
00643 cnl=lgnodes[j].cn[i+1][1];
00644 if (cnl==0){ ii++; }
00645 else{
00646 cm1=constrmat[ii][ns*2+0];
00647 cm2=constrmat[ii][ns*2+1];
00648 cn1=cn[ns][j][1];
00649 if (cn1>0) lv1=lv[cn1-1];
00650 else lv1=0.0;
00651 cn2=cn[ns][j][3];
00652 if (cn2>0) lv2=lv[cn2-1];
00653 else lv2=0.0;
00654 gv[cnl-1]+=cm1*lv1+cm2*lv2;
00655 ii++;
00656 }
00657
00658
00659 cnl=lgnodes[j].cn[i+1][2];
00660 if (cnl==0){ ii++; }
00661 else{
00662 cm1=constrmat[ii][ns*2+0];
00663 cn1=cn[ns][j][2];
00664 if (cn1>0) lv1=lv[cn1-1];
00665 else lv1=0.0;
00666 gv[cnl-1]+=cm1*lv1;
00667 ii++;
00668 }
00669
00670 cnl=lgnodes[j].cn[i+1][3];
00671 if (cnl==0){ ii++; }
00672 else{
00673 cm1=constrmat[ii][ns*2+0];
00674 cn1=cn[ns][j][5];
00675 if (cn1>0) lv1=lv[cn1-1];
00676 else lv1=0.0;
00677 gv[cnl-1]+=cm1*lv1;
00678 ii++;
00679 }
00680 }
00681 }
00682 }
00683
00684
00685 void lplate::locglob_auxiliary (double *gv,double *lv,long did,FILE *out)
00686 {
00687 long i,j;
00688 double s,*av;
00689
00690 av = new double [ngdof];
00691
00692
00693
00694
00695
00696
00697
00698
00699
00700
00701
00702
00703 nullvr (av,ngdof);
00704 j=0;
00705 for (i=adr[did];i<adr[did+1];i++){
00706 av[i]=lv[j]; j++;
00707 }
00708
00709
00710
00711
00712
00713
00714
00715
00716 for (i=0;i<nmult;i++){
00717 s=0.0;
00718 for (j=0;j<ngdof;j++){
00719 s+=dcm[i][j]*av[j];
00720 }
00721 gv[i]+=s;
00722 }
00723
00724 delete [] av;
00725 }
00726 void lplate::globloc_auxiliary (double *gv,double *lv,long did)
00727 {
00728 long i,j;
00729 double s,*av;
00730
00731 av = new double [ngdof];
00732
00733 nullvr (av,ngdof);
00734 for (i=0;i<ngdof;i++){
00735 s=0.0;
00736 for (j=0;j<nmult;j++){
00737 s+=dcm[j][i]*gv[j];
00738 }
00739 av[i]=s;
00740 }
00741
00742
00743
00744
00745
00746
00747
00748
00749
00750
00751
00752 j=0;
00753 for (i=adr[did];i<adr[did+1];i++){
00754 lv[j]=av[i]; j++;
00755 }
00756
00757 delete [] av;
00758 }
00759
00760
00761
00762
00763
00764
00765
00766
00767
00768
00769
00770
00771
00772
00773 void lplate::globloc (double *gv,double *lv,long ns)
00774 {
00775 long i,j,li,ii,cnl,cn1,cn2;
00776 double cm1,cm2,cgv;
00777
00778 li=ns-1;
00779 if (li<0) li=0;
00780 ii=li*nn*nnmult;
00781 for (i=li;i<nl-1;i++){
00782 for (j=0;j<nn;j++){
00783
00784 cnl=lgnodes[j].cn[i+1][0];
00785 if (cnl==0){ ii++; }
00786 else{
00787 cm1=constrmat[ii][ns*2+0];
00788 cm2=constrmat[ii][ns*2+1];
00789 cgv=gv[cnl-1];
00790
00791 cn1=cn[ns][j][0];
00792 if (cn1>0) lv[cn1-1]+=cgv*cm1;
00793 cn2=cn[ns][j][4];
00794 if (cn2>0) lv[cn2-1]+=cgv*cm2;
00795 ii++;
00796 }
00797
00798 cnl=lgnodes[j].cn[i+1][1];
00799 if (cnl==0){ ii++; }
00800 else{
00801 cm1=constrmat[ii][ns*2+0];
00802 cm2=constrmat[ii][ns*2+1];
00803 cgv=gv[cnl-1];
00804
00805 cn1=cn[ns][j][1];
00806 if (cn1>0) lv[cn1-1]+=cgv*cm1;
00807 cn2=cn[ns][j][3];
00808 if (cn2>0) lv[cn2-1]+=cgv*cm2;
00809 ii++;
00810 }
00811
00812 cnl=lgnodes[j].cn[i+1][2];
00813 if (cnl==0){ ii++; }
00814 else{
00815 cm1=constrmat[ii][ns*2+0];
00816 cgv=gv[cnl-1];
00817
00818 cn1=cn[ns][j][2];
00819 if (cn1>0) lv[cn1-1]+=cgv*cm1;
00820 ii++;
00821 }
00822
00823 cnl=lgnodes[j].cn[i+1][3];
00824 if (cnl==0){ ii++; }
00825 else{
00826 cm1=constrmat[ii][ns*2+0];
00827 cgv=gv[cnl-1];
00828
00829 cn1=cn[ns][j][5];
00830 if (cn1>0) lv[cn1-1]+=cgv*cm1;
00831 ii++;
00832 }
00833 }
00834 }
00835
00836 }
00837
00838
00839
00840
00841
00842
00843
00844
00845 void lplate::printconstrmat (FILE *out)
00846 {
00847 long i,j;
00848
00849 fprintf (out,"\n\n\n Constraint matrix \n");
00850 for (i=0;i<(nl-1)*nn*4;i++){
00851 fprintf (out,"\n row %4ld ",i);
00852 for (j=0;j<nl*2;j++){
00853 fprintf (out," %6.3e",constrmat[i][j]);
00854
00855 }
00856 }
00857 }
00858
00859
00860
00861
00862
00863
00864
00865
00866 void lplate::printconstrmat_dm (FILE *out)
00867 {
00868 long i,j;
00869
00870 fprintf (out,"\n\n\n Constraint matrix \n");
00871 for (i=0;i<nmult;i++){
00872 fprintf (out,"\n row %4ld ",i);
00873 for (j=0;j<nproc*ndof;j++){
00874 fprintf (out," %6.3e",dcm[i][j]);
00875
00876 }
00877 }
00878 }
00879
00880 void lplate::precondstep (double *dd,double *pp)
00881 {
00882 long i;
00883
00884 for (i=0;i<ndof;i++){
00885 pp[i]=dd[i];
00886 }
00887 }
00888
00889
00890
00891
00892
00893
00894
00895
00896
00897
00898
00899
00900
00901 void lplate::cg (gmatrix *gm,double *w,double *rhs,long *domproc,FILE *out)
00902 {
00903 long i,j;
00904 double nom,denom,alpha,beta,norrhs,pr1;
00905 double *d,*dd,*r,*p,*pp,*h;
00906 MPI_Status stat;
00907
00908 fprintf (stdout,"\n ndof %ld",ndof);
00909 fprintf (stdout,"\n maxndof %ld",maxndof);
00910 fprintf (stdout,"\n nmult %ld",nmult);
00911
00912 fprintf (out,"\n ndof %ld",ndof);
00913 fprintf (out,"\n maxndof %ld",maxndof);
00914 fprintf (out,"\n nmult %ld",nmult);
00915
00916
00917
00918 dd = new double [maxndof+1];
00919
00920 pp = new double [maxndof+1];
00921
00922 nom=0.0;
00923 for (i=0;i<ndof;i++){
00924 nom+=rhs[i]*rhs[i];
00925 }
00926
00927 if (myrank==0){
00928 r = new double [nmult];
00929 d = new double [nmult];
00930 p = new double [nmult];
00931 h = new double [nmult];
00932
00933
00934
00935 for (j=1;j<nproc;j++){
00936 nullvr (dd,maxndof+1);
00937 globloc (w,dd,domproc[j]);
00938
00939 MPI_Send (dd,maxndof,MPI_DOUBLE,j,ndom,MPI_COMM_WORLD);
00940 }
00941 nullvr (dd,maxndof+1);
00942 globloc (w,dd,domproc[0]);
00943
00944
00945
00946
00947 for (j=0;j<ndof;j++){
00948 dd[j]=rhs[j]-dd[j];
00949 }
00950
00951
00952 gm->back_substitution (pp,dd);
00953
00954
00955 nullvr (r,nmult);
00956 locglob (r,pp,domproc[0]);
00957
00958
00959 norrhs=nom;
00960
00961 for (j=1;j<nproc;j++){
00962 MPI_Recv(pp,maxndof+1,MPI_DOUBLE,MPI_ANY_SOURCE,MPI_ANY_TAG,MPI_COMM_WORLD,&stat);
00963 locglob (r,pp,domproc[stat.MPI_TAG]);
00964
00965 norrhs+=pp[maxndof];
00966 }
00967
00968 fprintf (out,"\n norm of right hand side vector %e,",norrhs);
00969 fprintf (stdout,"\n norm of right hand side vector %e,",norrhs);
00970
00971
00972 for (j=1;j<nproc;j++){
00973 nullvr (dd,maxndof+1);
00974 globloc (r,dd,domproc[j]);
00975 MPI_Send (dd,maxndof,MPI_DOUBLE,j,ndom,MPI_COMM_WORLD);
00976 }
00977 nullvr (dd,maxndof+1);
00978 globloc (r,dd,domproc[0]);
00979
00980 precondstep (dd,pp);
00981
00982 nullvr (h,nmult);
00983 locglob (h,pp,domproc[0]);
00984 for (j=1;j<nproc;j++){
00985 MPI_Recv(pp,maxndof,MPI_DOUBLE,MPI_ANY_SOURCE,MPI_ANY_TAG,MPI_COMM_WORLD,&stat);
00986 locglob (h,pp,domproc[stat.MPI_TAG]);
00987 }
00988
00989
00990 for (i=0;i<nmult;i++){
00991 d[i]=h[i];
00992 }
00993
00994
00995 nom = ss (r,h,nmult);
00996
00997
00998
00999
01000
01001 for (i=0;i<nicg;i++){
01002
01003
01004
01005
01006
01007
01008 for (j=1;j<nproc;j++){
01009 nullvr (dd,maxndof+1);
01010 globloc (d,dd,domproc[j]);
01011
01012 MPI_Send (dd,maxndof+1,MPI_DOUBLE,j,ndom,MPI_COMM_WORLD);
01013 }
01014
01015 nullvr (dd,maxndof+1);
01016 globloc (d,dd,domproc[0]);
01017
01018
01019
01020 gm->back_substitution (pp,dd);
01021
01022
01023 nullvr (p,nmult);
01024 locglob (p,pp,domproc[0]);
01025
01026
01027 for (j=1;j<nproc;j++){
01028 MPI_Recv(pp,maxndof+1,MPI_DOUBLE,MPI_ANY_SOURCE,MPI_ANY_TAG,MPI_COMM_WORLD,&stat);
01029 locglob (p,pp,domproc[stat.MPI_TAG]);
01030
01031 }
01032
01033
01034 denom = ss (d,p,nmult);
01035
01036 pr1=denom;
01037
01038 if (fabs(denom)<zero){
01039 fprintf (stderr,"\n\n zero denominator in modified conjugate gradient method (file %s, line %d).\n",__FILE__,__LINE__);
01040
01041 if (i!=nicg-1){
01042 dd[maxndof]=1.0;
01043 for (j=1;j<nproc;j++){
01044 MPI_Send (dd,maxndof+1,MPI_DOUBLE,j,0,MPI_COMM_WORLD);
01045 }
01046 }
01047 break;
01048 }
01049
01050
01051
01052
01053 alpha = nom/denom;
01054
01055
01056
01057
01058 for (j=0;j<nmult;j++){
01059 w[j]+=alpha*d[j];
01060 r[j]-=alpha*p[j];
01061 }
01062
01063 denom = nom;
01064 if (fabs(denom)<zero){
01065 fprintf (stderr,"\n\n zero denominator of beta in function (file %s, line %d).\n",__FILE__,__LINE__);
01066
01067 if (i!=nicg-1){
01068 dd[maxndof]=1.0;
01069 for (j=1;j<nproc;j++){
01070 MPI_Send (dd,maxndof+1,MPI_DOUBLE,j,myrank,MPI_COMM_WORLD);
01071 }
01072 }
01073 break;
01074 }
01075
01076
01077
01078 for (j=1;j<nproc;j++){
01079 nullvr (dd,maxndof+1);
01080 globloc (r,dd,domproc[j]);
01081 MPI_Send (dd,maxndof+1,MPI_DOUBLE,j,ndom,MPI_COMM_WORLD);
01082 }
01083 nullvr (dd,maxndof+1);
01084 globloc (r,dd,domproc[0]);
01085
01086 precondstep (dd,pp);
01087
01088 nullvr (h,nmult);
01089 locglob (h,pp,domproc[0]);
01090 for (j=1;j<nproc;j++){
01091 MPI_Recv(pp,maxndof+1,MPI_DOUBLE,MPI_ANY_SOURCE,MPI_ANY_TAG,MPI_COMM_WORLD,&stat);
01092 locglob (h,pp,domproc[stat.MPI_TAG]);
01093 }
01094
01095
01096 nom = ss (r,h,nmult);
01097
01098 fprintf (stdout,"\n iteration %4ld norm r %e mod. error %e alpha denom %e",i,sqrt(nom),sqrt(nom)/sqrt(norrhs),pr1);
01099 fprintf (out,"\n iteration %4ld norm r %e mod. error %e alpha denom %e",i,sqrt(nom),sqrt(nom)/sqrt(norrhs),pr1);
01100
01101
01102 if (sqrt(nom)/sqrt(norrhs)<errcg){
01103 if (i!=nicg-1){
01104 dd[maxndof]=1.0;
01105 for (j=1;j<nproc;j++){
01106 MPI_Send (dd,maxndof+1,MPI_DOUBLE,j,myrank,MPI_COMM_WORLD);
01107 }
01108 }
01109 break;
01110 }
01111
01112
01113
01114 beta = nom/denom;
01115
01116
01117 for (j=0;j<nmult;j++){
01118 d[j]=beta*d[j]+h[j];
01119 }
01120
01121 }
01122
01123
01124 anicg=i; aerrcg=nom;
01125
01126 fprintf (out,"\n\n\n\n kontrola Lagrangeovych multiplikatoru \n");
01127 for (i=0;i<nmult;i++){
01128 fprintf (out,"\n lagr. mult %4ld %e",i,w[i]);
01129 }
01130
01131 }
01132 else{
01133
01134 MPI_Recv (dd,maxndof,MPI_DOUBLE,MPI_ANY_SOURCE,MPI_ANY_TAG,MPI_COMM_WORLD,&stat);
01135
01136
01137 for (j=0;j<ndof;j++){
01138 dd[j]=rhs[j]-dd[j];
01139 }
01140
01141
01142 gm->back_substitution (pp,dd);
01143
01144
01145 pp[maxndof]=nom;
01146
01147 MPI_Send (pp,maxndof+1,MPI_DOUBLE,0,ndom,MPI_COMM_WORLD);
01148
01149
01150
01151 MPI_Recv (dd,maxndof,MPI_DOUBLE,MPI_ANY_SOURCE,MPI_ANY_TAG,MPI_COMM_WORLD,&stat);
01152 precondstep (dd,pp);
01153 MPI_Send (pp,maxndof,MPI_DOUBLE,0,ndom,MPI_COMM_WORLD);
01154
01155 for (i=0;i<nicg;i++){
01156
01157 MPI_Recv (dd,maxndof+1,MPI_DOUBLE,MPI_ANY_SOURCE,MPI_ANY_TAG,MPI_COMM_WORLD,&stat);
01158 if (dd[maxndof]>0.5){ break; }
01159
01160
01161 gm->back_substitution (pp,dd);
01162
01163
01164 MPI_Send (pp,maxndof+1,MPI_DOUBLE,0,myrank,MPI_COMM_WORLD);
01165
01166
01167 MPI_Recv (dd,maxndof+1,MPI_DOUBLE,MPI_ANY_SOURCE,MPI_ANY_TAG,MPI_COMM_WORLD,&stat);
01168 if (dd[maxndof]>0.5){ break; }
01169 precondstep (dd,pp);
01170 MPI_Send (pp,maxndof+1,MPI_DOUBLE,0,ndom,MPI_COMM_WORLD);
01171
01172 }
01173
01174 }
01175
01176 if (myrank==0){
01177 delete [] r;
01178 delete [] d;
01179 delete [] p;
01180 }
01181
01182
01183
01184
01185 delete [] dd;
01186 delete [] pp;
01187
01188 }
01189
01190
01191
01192
01193
01194
01195
01196
01197
01198
01199
01200
01201 void lplate::nodaldisplacements (gmatrix *gm,double *lhs,double *rhs,double *w,long *domproc)
01202 {
01203 long j;
01204 double *dd;
01205 MPI_Status stat;
01206
01207
01208 dd = new double [maxndof+1];
01209
01210 if (myrank==0){
01211
01212 for (j=1;j<nproc;j++){
01213 nullvr (dd,maxndof+1);
01214 globloc (w,dd,domproc[j]);
01215
01216 MPI_Send (dd,maxndof,MPI_DOUBLE,j,ndom,MPI_COMM_WORLD);
01217 }
01218
01219 nullvr (dd,maxndof+1);
01220 globloc (w,dd,domproc[0]);
01221
01222
01223
01224 for (j=0;j<ndof;j++){
01225 dd[j]=rhs[j]-dd[j];
01226 }
01227
01228
01229 gm->back_substitution (lhs,dd);
01230
01231 }
01232 else{
01233
01234 MPI_Recv (dd,maxndof,MPI_DOUBLE,MPI_ANY_SOURCE,MPI_ANY_TAG,MPI_COMM_WORLD,&stat);
01235
01236
01237 for (j=0;j<ndof;j++){
01238 dd[j]=rhs[j]-dd[j];
01239 }
01240
01241
01242 gm->back_substitution (lhs,dd);
01243
01244 }
01245
01246 delete [] dd;
01247 }
01248
01249
01250
01251
01252
01253
01254
01255
01256
01257
01258 void lplate::solve_system (gtopology *top,gmatrix *gm,
01259 long *domproc,double *lhs,double *rhs,FILE *out)
01260 {
01261 double *w;
01262
01263
01264
01265
01266
01267
01268
01269
01270
01271
01272 if (myrank==0){
01273 w = new double [nmult];
01274 nullvr (w,nmult);
01275 }
01276
01277
01278 gm->decompose_matrix ();
01279
01280
01281
01282
01283
01284 if (myrank==0) orthonormalization ();
01285
01286
01287
01288
01289
01290
01291 cg (gm,w,rhs,domproc,out);
01292
01293
01294 nodaldisplacements (gm,lhs,rhs,w,domproc);
01295
01296 if (myrank==0)
01297 delete [] w;
01298
01299 }