00001 #include "lplate.h"
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012 lplate::lplate (int np,int mr,long 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 nullv (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 nullv (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 (gmatrix *gm,double *dd,double *pp)
00881 {
00882
00883
00884
00885
00886
00887
00888
00889
00890 gm -> decompgmxv (dd,pp);
00891 }
00892
00893
00894
00895
00896
00897
00898
00899
00900
00901
00902
00903
00904
00905 void lplate::cg (gmatrix *gm,double *w,double *rhs,long *domproc,FILE *out)
00906 {
00907 long i,j;
00908 double nom,denom,alpha,beta,norrhs,pr1;
00909 double *d,*dd,*r,*p,*pp;
00910 MPI_Status stat;
00911
00912 fprintf (stdout,"\n ndof %ld",ndof);
00913 fprintf (stdout,"\n maxndof %ld",maxndof);
00914 fprintf (stdout,"\n nmult %ld",nmult);
00915
00916 fprintf (out,"\n ndof %ld",ndof);
00917 fprintf (out,"\n maxndof %ld",maxndof);
00918 fprintf (out,"\n nmult %ld",nmult);
00919
00920
00921
00922 dd = new double [maxndof+1];
00923
00924 pp = new double [maxndof+1];
00925
00926 nom=0.0;
00927 for (i=0;i<ndof;i++){
00928 nom+=rhs[i]*rhs[i];
00929 }
00930
00931 if (myrank==0){
00932 r = new double [nmult];
00933 d = new double [nmult];
00934 p = new double [nmult];
00935
00936 nullv (w,nmult);
00937
00938
00939 for (j=1;j<nproc;j++){
00940 nullv (dd,maxndof+1);
00941 globloc (w,dd,domproc[j]);
00942
00943 MPI_Send (dd,maxndof,MPI_DOUBLE,j,ndom,MPI_COMM_WORLD);
00944 }
00945 nullv (dd,maxndof+1);
00946 globloc (w,dd,domproc[0]);
00947
00948
00949
00950 for (j=0;j<ndof;j++){
00951 dd[j]=rhs[j]-dd[j];
00952 }
00953
00954
00955 gm->back_substitution (pp,dd);
00956
00957
00958 nullv (r,nmult);
00959 locglob (r,pp,domproc[0]);
00960
00961
00962 norrhs=nom;
00963
00964 for (j=1;j<nproc;j++){
00965 MPI_Recv(pp,maxndof+1,MPI_DOUBLE,MPI_ANY_SOURCE,MPI_ANY_TAG,MPI_COMM_WORLD,&stat);
00966 locglob (r,pp,domproc[stat.MPI_TAG]);
00967
00968 norrhs+=pp[maxndof];
00969 }
00970
00971 fprintf (out,"\n norm of right hand side vector %e,",norrhs);
00972 fprintf (stdout,"\n norm of right hand side vector %e,",norrhs);
00973
00974
00975 for (i=0;i<nmult;i++){
00976 d[i]=r[i];
00977 }
00978
00979
00980 nom = ss (r,r,nmult);
00981
00982
00983
00984
00985
00986 for (i=0;i<nicg;i++){
00987
00988
00989
00990
00991
00992
00993 for (j=1;j<nproc;j++){
00994 nullv (dd,maxndof+1);
00995 globloc (d,dd,domproc[j]);
00996
00997 MPI_Send (dd,maxndof+1,MPI_DOUBLE,j,ndom,MPI_COMM_WORLD);
00998 }
00999
01000 nullv (dd,maxndof+1);
01001 globloc (d,dd,domproc[0]);
01002
01003
01004
01005 gm->back_substitution (pp,dd);
01006
01007
01008
01009 nullv (p,nmult);
01010 locglob (p,pp,domproc[0]);
01011
01012
01013 for (j=1;j<nproc;j++){
01014 MPI_Recv(pp,maxndof+1,MPI_DOUBLE,MPI_ANY_SOURCE,MPI_ANY_TAG,MPI_COMM_WORLD,&stat);
01015 locglob (p,pp,domproc[stat.MPI_TAG]);
01016
01017 }
01018
01019
01020 denom = ss (d,p,nmult);
01021
01022 pr1=denom;
01023
01024 if (fabs(denom)<zero){
01025 fprintf (stderr,"\n\n zero denominator in modified conjugate gradient method (file %s, line %d).\n",__FILE__,__LINE__);
01026
01027 if (i!=nicg-1){
01028 dd[maxndof]=1.0;
01029 for (j=1;j<nproc;j++){
01030 MPI_Send (dd,maxndof+1,MPI_DOUBLE,j,0,MPI_COMM_WORLD);
01031 }
01032 }
01033 break;
01034 }
01035
01036
01037
01038
01039 alpha = nom/denom;
01040
01041
01042
01043
01044 for (j=0;j<nmult;j++){
01045 w[j]+=alpha*d[j];
01046 r[j]-=alpha*p[j];
01047 }
01048
01049 denom = nom;
01050 if (fabs(denom)<zero){
01051 fprintf (stderr,"\n\n zero denominator of beta in function (file %s, line %d).\n",__FILE__,__LINE__);
01052
01053 if (i!=nicg-1){
01054 dd[maxndof]=1.0;
01055 for (j=1;j<nproc;j++){
01056 MPI_Send (dd,maxndof+1,MPI_DOUBLE,j,myrank,MPI_COMM_WORLD);
01057 }
01058 }
01059 break;
01060 }
01061
01062 nom = ss (r,r,nmult);
01063
01064 fprintf (stdout,"\n iteration %4ld norm r %e mod. error %e alpha denom %e",i,sqrt(nom),sqrt(nom)/sqrt(norrhs),pr1);
01065 fprintf (out,"\n iteration %4ld norm r %e mod. error %e alpha denom %e",i,sqrt(nom),sqrt(nom)/sqrt(norrhs),pr1);
01066
01067
01068 if (sqrt(nom)/sqrt(norrhs)<errcg){
01069 if (i!=nicg-1){
01070 dd[maxndof]=1.0;
01071 for (j=1;j<nproc;j++){
01072 MPI_Send (dd,maxndof+1,MPI_DOUBLE,j,myrank,MPI_COMM_WORLD);
01073 }
01074 }
01075 break;
01076 }
01077
01078
01079
01080 beta = nom/denom;
01081
01082
01083 for (j=0;j<nmult;j++){
01084 d[j]=beta*d[j]+r[j];
01085 }
01086
01087 }
01088
01089
01090 anicg=i; aerrcg=nom;
01091
01092 fprintf (out,"\n\n\n\n kontrola Lagrangeovych multiplikatoru \n");
01093 for (i=0;i<nmult;i++){
01094 fprintf (out,"\n lagr. mult %4ld %e",i,w[i]);
01095 }
01096
01097 }
01098 else{
01099
01100 MPI_Recv (dd,maxndof,MPI_DOUBLE,MPI_ANY_SOURCE,MPI_ANY_TAG,MPI_COMM_WORLD,&stat);
01101
01102
01103 for (j=0;j<ndof;j++){
01104 dd[j]=rhs[j]-dd[j];
01105 }
01106
01107
01108
01109 gm->back_substitution (pp,dd);
01110
01111
01112 pp[maxndof]=nom;
01113
01114 MPI_Send (pp,maxndof+1,MPI_DOUBLE,0,ndom,MPI_COMM_WORLD);
01115
01116
01117 for (i=0;i<nicg;i++){
01118
01119
01120 MPI_Recv (dd,maxndof+1,MPI_DOUBLE,MPI_ANY_SOURCE,MPI_ANY_TAG,MPI_COMM_WORLD,&stat);
01121 if (dd[maxndof]>0.5){ break; }
01122
01123
01124 gm->back_substitution (pp,dd);
01125
01126
01127 MPI_Send (pp,maxndof+1,MPI_DOUBLE,0,myrank,MPI_COMM_WORLD);
01128
01129 }
01130
01131 }
01132
01133 if (myrank==0){
01134 delete [] r;
01135 delete [] d;
01136 delete [] p;
01137 }
01138
01139
01140
01141
01142 delete [] dd;
01143 delete [] pp;
01144
01145 }
01146
01147
01148
01149
01150
01151
01152
01153
01154
01155
01156
01157
01158 void lplate::pcg (gmatrix *gm,double *w,double *rhs,long *domproc,FILE *out)
01159 {
01160 long i,j;
01161 double nom,denom,alpha,beta,norrhs,pr1;
01162 double *d,*dd,*r,*p,*pp,*h;
01163 MPI_Status stat;
01164
01165 fprintf (stdout,"\n ndof %ld",ndof);
01166 fprintf (stdout,"\n maxndof %ld",maxndof);
01167 fprintf (stdout,"\n nmult %ld",nmult);
01168
01169 fprintf (out,"\n ndof %ld",ndof);
01170 fprintf (out,"\n maxndof %ld",maxndof);
01171 fprintf (out,"\n nmult %ld",nmult);
01172
01173
01174
01175 dd = new double [maxndof+1];
01176
01177 pp = new double [maxndof+1];
01178
01179 nom=0.0;
01180 for (i=0;i<ndof;i++){
01181 nom+=rhs[i]*rhs[i];
01182 }
01183
01184 if (myrank==0){
01185 r = new double [nmult];
01186 d = new double [nmult];
01187 p = new double [nmult];
01188 h = new double [nmult];
01189
01190
01191
01192 for (j=1;j<nproc;j++){
01193 nullv (dd,maxndof+1);
01194 globloc (w,dd,domproc[j]);
01195
01196 MPI_Send (dd,maxndof,MPI_DOUBLE,j,ndom,MPI_COMM_WORLD);
01197 }
01198 nullv (dd,maxndof+1);
01199 globloc (w,dd,domproc[0]);
01200
01201
01202
01203
01204 for (j=0;j<ndof;j++){
01205 dd[j]=rhs[j]-dd[j];
01206 }
01207
01208
01209 gm->back_substitution (pp,dd);
01210
01211
01212 nullv (r,nmult);
01213 locglob (r,pp,domproc[0]);
01214
01215
01216 norrhs=nom;
01217
01218 for (j=1;j<nproc;j++){
01219 MPI_Recv(pp,maxndof+1,MPI_DOUBLE,MPI_ANY_SOURCE,MPI_ANY_TAG,MPI_COMM_WORLD,&stat);
01220 locglob (r,pp,domproc[stat.MPI_TAG]);
01221
01222 norrhs+=pp[maxndof];
01223 }
01224
01225 fprintf (out,"\n norm of right hand side vector %e,",norrhs);
01226 fprintf (stdout,"\n norm of right hand side vector %e,",norrhs);
01227
01228
01229 for (j=1;j<nproc;j++){
01230 nullv (dd,maxndof+1);
01231 globloc (r,dd,domproc[j]);
01232 MPI_Send (dd,maxndof,MPI_DOUBLE,j,ndom,MPI_COMM_WORLD);
01233 }
01234 nullv (dd,maxndof+1);
01235 globloc (r,dd,domproc[0]);
01236
01237 precondstep (gm,dd,pp);
01238
01239 nullv (h,nmult);
01240 locglob (h,pp,domproc[0]);
01241 for (j=1;j<nproc;j++){
01242 MPI_Recv(pp,maxndof,MPI_DOUBLE,MPI_ANY_SOURCE,MPI_ANY_TAG,MPI_COMM_WORLD,&stat);
01243 locglob (h,pp,domproc[stat.MPI_TAG]);
01244 }
01245
01246
01247 for (i=0;i<nmult;i++){
01248 d[i]=h[i];
01249 }
01250
01251
01252 nom = ss (r,h,nmult);
01253
01254
01255
01256
01257
01258 for (i=0;i<nicg;i++){
01259
01260
01261
01262
01263
01264
01265 for (j=1;j<nproc;j++){
01266 nullv (dd,maxndof+1);
01267 globloc (d,dd,domproc[j]);
01268
01269 MPI_Send (dd,maxndof+1,MPI_DOUBLE,j,ndom,MPI_COMM_WORLD);
01270 }
01271
01272 nullv (dd,maxndof+1);
01273 globloc (d,dd,domproc[0]);
01274
01275
01276
01277 gm->back_substitution (pp,dd);
01278
01279
01280 nullv (p,nmult);
01281 locglob (p,pp,domproc[0]);
01282
01283
01284 for (j=1;j<nproc;j++){
01285 MPI_Recv(pp,maxndof+1,MPI_DOUBLE,MPI_ANY_SOURCE,MPI_ANY_TAG,MPI_COMM_WORLD,&stat);
01286 locglob (p,pp,domproc[stat.MPI_TAG]);
01287
01288 }
01289
01290
01291 denom = ss (d,p,nmult);
01292
01293 pr1=denom;
01294
01295 if (fabs(denom)<zero){
01296 fprintf (stderr,"\n\n zero denominator in modified conjugate gradient method (file %s, line %d).\n",__FILE__,__LINE__);
01297
01298 if (i!=nicg-1){
01299 dd[maxndof]=1.0;
01300 for (j=1;j<nproc;j++){
01301 MPI_Send (dd,maxndof+1,MPI_DOUBLE,j,0,MPI_COMM_WORLD);
01302 }
01303 }
01304 break;
01305 }
01306
01307
01308
01309
01310 alpha = nom/denom;
01311
01312
01313
01314
01315 for (j=0;j<nmult;j++){
01316 w[j]+=alpha*d[j];
01317 r[j]-=alpha*p[j];
01318 }
01319
01320 denom = nom;
01321 if (fabs(denom)<zero){
01322 fprintf (stderr,"\n\n zero denominator of beta in function (file %s, line %d).\n",__FILE__,__LINE__);
01323
01324 if (i!=nicg-1){
01325 dd[maxndof]=1.0;
01326 for (j=1;j<nproc;j++){
01327 MPI_Send (dd,maxndof+1,MPI_DOUBLE,j,myrank,MPI_COMM_WORLD);
01328 }
01329 }
01330 break;
01331 }
01332
01333
01334
01335 for (j=1;j<nproc;j++){
01336 nullv (dd,maxndof+1);
01337 globloc (r,dd,domproc[j]);
01338 MPI_Send (dd,maxndof+1,MPI_DOUBLE,j,ndom,MPI_COMM_WORLD);
01339 }
01340 nullv (dd,maxndof+1);
01341 globloc (r,dd,domproc[0]);
01342
01343 precondstep (gm,dd,pp);
01344
01345 nullv (h,nmult);
01346 locglob (h,pp,domproc[0]);
01347 for (j=1;j<nproc;j++){
01348 MPI_Recv(pp,maxndof+1,MPI_DOUBLE,MPI_ANY_SOURCE,MPI_ANY_TAG,MPI_COMM_WORLD,&stat);
01349 locglob (h,pp,domproc[stat.MPI_TAG]);
01350 }
01351
01352
01353 nom = ss (r,h,nmult);
01354
01355 fprintf (stdout,"\n iteration %4ld norm r %e mod. error %e alpha denom %e",i,sqrt(nom),sqrt(nom)/sqrt(norrhs),pr1);
01356 fprintf (out,"\n iteration %4ld norm r %e mod. error %e alpha denom %e",i,sqrt(nom),sqrt(nom)/sqrt(norrhs),pr1);
01357
01358
01359 if (sqrt(nom)/sqrt(norrhs)<errcg){
01360 if (i!=nicg-1){
01361 dd[maxndof]=1.0;
01362 for (j=1;j<nproc;j++){
01363 MPI_Send (dd,maxndof+1,MPI_DOUBLE,j,myrank,MPI_COMM_WORLD);
01364 }
01365 }
01366 break;
01367 }
01368
01369
01370
01371 beta = nom/denom;
01372
01373
01374 for (j=0;j<nmult;j++){
01375 d[j]=beta*d[j]+h[j];
01376 }
01377
01378 }
01379
01380
01381 anicg=i; aerrcg=nom;
01382
01383 fprintf (out,"\n\n\n\n kontrola Lagrangeovych multiplikatoru \n");
01384 for (i=0;i<nmult;i++){
01385 fprintf (out,"\n lagr. mult %4ld %e",i,w[i]);
01386 }
01387
01388 }
01389 else{
01390
01391 MPI_Recv (dd,maxndof,MPI_DOUBLE,MPI_ANY_SOURCE,MPI_ANY_TAG,MPI_COMM_WORLD,&stat);
01392
01393
01394 for (j=0;j<ndof;j++){
01395 dd[j]=rhs[j]-dd[j];
01396 }
01397
01398
01399 gm->back_substitution (pp,dd);
01400
01401
01402 pp[maxndof]=nom;
01403
01404 MPI_Send (pp,maxndof+1,MPI_DOUBLE,0,ndom,MPI_COMM_WORLD);
01405
01406
01407
01408 MPI_Recv (dd,maxndof,MPI_DOUBLE,MPI_ANY_SOURCE,MPI_ANY_TAG,MPI_COMM_WORLD,&stat);
01409 precondstep (gm,dd,pp);
01410 MPI_Send (pp,maxndof,MPI_DOUBLE,0,ndom,MPI_COMM_WORLD);
01411
01412 for (i=0;i<nicg;i++){
01413
01414 MPI_Recv (dd,maxndof+1,MPI_DOUBLE,MPI_ANY_SOURCE,MPI_ANY_TAG,MPI_COMM_WORLD,&stat);
01415 if (dd[maxndof]>0.5){ break; }
01416
01417
01418 gm->back_substitution (pp,dd);
01419
01420
01421 MPI_Send (pp,maxndof+1,MPI_DOUBLE,0,myrank,MPI_COMM_WORLD);
01422
01423
01424 MPI_Recv (dd,maxndof+1,MPI_DOUBLE,MPI_ANY_SOURCE,MPI_ANY_TAG,MPI_COMM_WORLD,&stat);
01425 if (dd[maxndof]>0.5){ break; }
01426 precondstep (gm,dd,pp);
01427 MPI_Send (pp,maxndof+1,MPI_DOUBLE,0,ndom,MPI_COMM_WORLD);
01428
01429 }
01430
01431 }
01432
01433 if (myrank==0){
01434 delete [] r;
01435 delete [] d;
01436 delete [] p;
01437 }
01438
01439
01440
01441
01442 delete [] dd;
01443 delete [] pp;
01444
01445 }
01446
01447
01448
01449
01450
01451
01452
01453
01454
01455
01456
01457
01458 void lplate::nodaldisplacements (gmatrix *gm,double *lhs,double *rhs,double *w,long *domproc)
01459 {
01460 long j;
01461 double *dd;
01462 MPI_Status stat;
01463
01464
01465 dd = new double [maxndof+1];
01466
01467 if (myrank==0){
01468
01469 for (j=1;j<nproc;j++){
01470 nullv (dd,maxndof+1);
01471 globloc (w,dd,domproc[j]);
01472
01473 MPI_Send (dd,maxndof,MPI_DOUBLE,j,ndom,MPI_COMM_WORLD);
01474 }
01475
01476 nullv (dd,maxndof+1);
01477 globloc (w,dd,domproc[0]);
01478
01479
01480
01481 for (j=0;j<ndof;j++){
01482 dd[j]=rhs[j]-dd[j];
01483 }
01484
01485
01486 gm->back_substitution (lhs,dd);
01487
01488 }
01489 else{
01490
01491 MPI_Recv (dd,maxndof,MPI_DOUBLE,MPI_ANY_SOURCE,MPI_ANY_TAG,MPI_COMM_WORLD,&stat);
01492
01493
01494 for (j=0;j<ndof;j++){
01495 dd[j]=rhs[j]-dd[j];
01496 }
01497
01498
01499 gm->back_substitution (lhs,dd);
01500
01501 }
01502
01503 delete [] dd;
01504 }
01505
01506
01507
01508
01509
01510
01511
01512
01513
01514
01515 void lplate::solve_system (gtopology *top,gmatrix *gm,
01516 long *domproc,double *lhs,double *rhs,FILE *out)
01517 {
01518 double *w;
01519 time_t t1,t2,t3,t4,t5;
01520
01521
01522
01523
01524
01525
01526
01527
01528
01529
01530 if (myrank==0){
01531 w = new double [nmult];
01532 nullv (w,nmult);
01533 }
01534
01535 t1 = time (NULL);
01536
01537
01538 gm->decompose_matrix ();
01539
01540 t2 = time (NULL);
01541
01542
01543
01544
01545
01546 if (myrank==0) orthonormalization ();
01547
01548
01549
01550
01551
01552 t3 = time (NULL);
01553
01554
01555 cg (gm,w,rhs,domproc,out);
01556
01557
01558 t4 = time (NULL);
01559
01560
01561 nodaldisplacements (gm,lhs,rhs,w,domproc);
01562
01563 t5 = time (NULL);
01564
01565 MPI_Status stat;
01566 long i,j,*buff;
01567 buff = new long[4];
01568 buff[0]=t2-t1;
01569 buff[1]=t3-t2;
01570 buff[2]=t4-t3;
01571 buff[3]=t5-t4;
01572 if (myrank==0){
01573
01574 fprintf (out,"\n\n\n\n\n");
01575 j=domproc[myrank];
01576 fprintf (out,"\n\n\n Domain %ld",j);
01577 fprintf (out,"\n time of condensation %ld",buff[0]);
01578 fprintf (out,"\n time of orthonormalization %ld",buff[1]);
01579 fprintf (out,"\n time of CG method %ld",buff[2]);
01580 fprintf (out,"\n time of displacement computation %ld",buff[3]);
01581
01582 j=domproc[myrank];
01583 fprintf (stdout,"\n\n\n Domain %ld",j);
01584 fprintf (out,"\n time of condensation %ld",buff[0]);
01585 fprintf (out,"\n time of orthonormalization %ld",buff[1]);
01586 fprintf (out,"\n time of CG method %ld",buff[2]);
01587 fprintf (out,"\n time of displacement computation %ld",buff[3]);
01588
01589
01590 for (i=1;i<nproc;i++){
01591 MPI_Recv (buff,4,MPI_LONG,MPI_ANY_SOURCE,MPI_ANY_TAG,MPI_COMM_WORLD,&stat);
01592 j=domproc[stat.MPI_TAG];
01593 fprintf (out,"\n\n\n Domain %ld",j);
01594 fprintf (out,"\n time of condensation %ld",buff[0]);
01595 fprintf (out,"\n time of orthonormalization %ld",buff[1]);
01596 fprintf (out,"\n time of CG method %ld",buff[2]);
01597 fprintf (out,"\n time of displacement computation %ld",buff[3]);
01598
01599 fprintf (stdout,"\n\n\n Domain %ld",j);
01600 fprintf (out,"\n time of condensation %ld",buff[0]);
01601 fprintf (out,"\n time of orthonormalization %ld",buff[1]);
01602 fprintf (out,"\n time of CG method %ld",buff[2]);
01603 fprintf (out,"\n time of displacement computation %ld",buff[3]);
01604 }
01605 }
01606 else{
01607 MPI_Send(buff,4,MPI_LONG,0,myrank,MPI_COMM_WORLD);
01608 }
01609 MPI_Barrier (MPI_COMM_WORLD);
01610 delete [] buff;
01611
01612 if (myrank==0)
01613 delete [] w;
01614
01615 }