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