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*4];
00054 for (i=0;i<(nl-1)*nn*4;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 for (i=0;i<nl;i++){
00070 delete [] thick[i];
00071 }
00072 delete [] thick;
00073
00074 for (i=0;i<nl;i++){
00075 for (j=0;j<nn;j++){
00076 delete [] cn[i][j];
00077 }
00078 delete [] cn[i];
00079 }
00080 delete [] cn;
00081
00082 for (i=0;i<(nl-1)*nn*4;i++){
00083 delete [] constrmat[i];
00084 }
00085 delete [] constrmat;
00086
00087 }
00088
00089
00090
00091
00092
00093
00094
00095
00096
00097
00098
00099
00100
00101
00102
00103
00104 void lplate::globcnnum_lpp (gtopology *top,long *domproc,FILE *out)
00105 {
00106 long i,j,k,buffsize;
00107 long *buff;
00108 MPI_Status stat;
00109
00110
00111 buffsize = nn*ndofn;
00112 buff = new long [buffsize];
00113
00114 k=0;
00115 for (i=0;i<nn;i++){
00116 for (j=0;j<ndofn;j++){
00117 buff[k]=top->gnodes[i].cn[j];
00118 k++;
00119 }
00120 }
00121
00122 if (myrank==0){
00123 assemble_cn (buff,ndom);
00124
00125 for (i=1;i<nproc;i++){
00126 MPI_Recv (buff,buffsize,MPI_LONG,MPI_ANY_SOURCE,MPI_ANY_TAG,MPI_COMM_WORLD,&stat);
00127 assemble_cn (buff,domproc[stat.MPI_TAG]);
00128 }
00129
00130
00131 maxndofdom ();
00132
00133 for (i=1;i<nproc;i++){
00134 MPI_Send (&maxndof,1,MPI_LONG,i,myrank,MPI_COMM_WORLD);
00135 }
00136
00137
00138 fprintf (out,"\n\n\n kontrola sebranych kodovych cisel \n");
00139 for (i=0;i<nproc;i++){
00140 fprintf (out,"\n\n DOMAIN number %ld",i);
00141 for (j=0;j<nn;j++){
00142 fprintf (out,"\n node %4ld ",j);
00143 for (k=0;k<ndofn;k++){
00144 fprintf (out," %ld",cn[i][j][k]);
00145 }
00146 }
00147 }
00148
00149
00150
00151 codenum_mult (nmult);
00152
00153
00154 fprintf (out,"\n\n\n kontrola nagenerovanych kodovych cisel \n");
00155 for (i=0;i<nn;i++){
00156 fprintf (out,"\n layered node number %4ld ",i);
00157 for (j=0;j<=nl;j++){
00158 fprintf (out," |%ld| ",j);
00159 for (k=0;k<nnmult;k++){
00160 fprintf (out," %ld",lgnodes[i].cn[j][k]);
00161 }
00162 }
00163 }
00164
00165 }
00166 else{
00167 MPI_Send (buff,buffsize,MPI_LONG,0,myrank,MPI_COMM_WORLD);
00168 MPI_Recv (&maxndof,1,MPI_LONG,MPI_ANY_SOURCE,MPI_ANY_TAG,MPI_COMM_WORLD,&stat);
00169 }
00170 MPI_Barrier (MPI_COMM_WORLD);
00171 delete [] buff;
00172
00173 }
00174
00175
00176
00177
00178
00179
00180
00181
00182
00183
00184
00185
00186 void lplate::constraintmat (double *th,long *domproc,FILE *out)
00187 {
00188 long i;
00189 double *rbuff;
00190 MPI_Status stat;
00191
00192
00193 rbuff = new double [nn];
00194 for (i=0;i<nn;i++){
00195 rbuff[i]=th[i];
00196 }
00197
00198 if (myrank==0){
00199 assemble_thicknesses (rbuff,ndom);
00200
00201 for (i=1;i<nproc;i++){
00202 MPI_Recv (rbuff,nn,MPI_DOUBLE,MPI_ANY_SOURCE,MPI_ANY_TAG,MPI_COMM_WORLD,&stat);
00203 assemble_thicknesses (rbuff,domproc[stat.MPI_TAG]);
00204 }
00205 }
00206 else{
00207 MPI_Send (rbuff,nn,MPI_DOUBLE,0,myrank,MPI_COMM_WORLD);
00208 }
00209 MPI_Barrier (MPI_COMM_WORLD);
00210 delete [] rbuff;
00211
00212
00213
00214 if (myrank==0){
00215 assemble_constr ();
00216
00217
00218
00219 }
00220 }
00221
00222
00223
00224
00225
00226
00227
00228 void lplate::maxndofdom ()
00229 {
00230 long i,j,k;
00231
00232 maxndof=0;
00233 for (i=0;i<nproc;i++){
00234 for (j=0;j<nn;j++){
00235 for (k=0;k<ndofn;k++){
00236 if (cn[i][j][k]>maxndof) maxndof=cn[i][j][k];
00237 }
00238 }
00239 }
00240 }
00241
00242
00243
00244
00245
00246
00247
00248
00249
00250
00251 void lplate::assemble_cn (long *buff,long ndom)
00252 {
00253 long i,j,k;
00254
00255 k=0;
00256 for (i=0;i<nn;i++){
00257 for (j=0;j<ndofn;j++){
00258 cn[ndom][i][j]=buff[k];
00259 k++;
00260 }
00261 }
00262 }
00263
00264
00265
00266
00267
00268
00269
00270
00271
00272
00273 void lplate::assemble_thicknesses (double *buff,long ndom)
00274 {
00275 long i;
00276
00277 for (i=0;i<nn;i++){
00278 thick[ndom][i]=buff[i];
00279 }
00280 }
00281
00282
00283
00284
00285
00286
00287
00288
00289
00290 void lplate::codenum_mult (long &nm)
00291 {
00292 long i,j;
00293
00294 nm=1;
00295 for (i=1;i<nl;i++){
00296 for (j=0;j<nn;j++){
00297 if (cn[i-1][j][0]>0 || cn[i-1][j][4]>0 || cn[i][j][0]>0 || cn[i][j][4]>0){
00298 lgnodes[j].cn[i][0]=nm; nm++;
00299 }
00300 if (cn[i-1][j][1]>0 || cn[i-1][j][3]>0 || cn[i][j][1]>0 || cn[i][j][3]>0){
00301 lgnodes[j].cn[i][1]=nm; nm++;
00302 }
00303 if (cn[i-1][j][2]>0 || cn[i][j][2]>0){
00304 lgnodes[j].cn[i][2]=nm; nm++;
00305 }
00306 if (cn[i-1][j][5]>0 || cn[i][j][5]>0){
00307 lgnodes[j].cn[i][3]=nm; nm++;
00308 }
00309 }
00310 }
00311 nm--;
00312 }
00313
00314
00315
00316
00317
00318
00319
00320
00321
00322
00323
00324
00325
00326
00327
00328 void lplate::assemble_constr ()
00329 {
00330 long i,j,ii;
00331 double lthickness,uthickness;
00332
00333 ii=0;
00334 for (i=1;i<nl;i++){
00335 for (j=0;j<nn;j++){
00336 lthickness = thick[i-1][j];
00337 uthickness = thick[i][j];
00338
00339 if (cn[i-1][j][0]>0) constrmat[ii][(i-1)*2+0]=-1.0;
00340 if (cn[i-1][j][4]>0) constrmat[ii][(i-1)*2+1]=lthickness/(-2.0);
00341 if (cn[i][j][0]>0) constrmat[ii][i*2+0]=1.0;
00342 if (cn[i][j][4]>0) constrmat[ii][i*2+1]=uthickness/(-2.0);
00343 ii++;
00344
00345 if (cn[i-1][j][1]>0) constrmat[ii][(i-1)*2+0]=-1.0;
00346 if (cn[i-1][j][3]>0) constrmat[ii][(i-1)*2+1]=lthickness/2.0;
00347 if (cn[i][j][1]>0) constrmat[ii][i*2+0]=1.0;
00348 if (cn[i][j][3]>0) constrmat[ii][i*2+1]=uthickness/2.0;
00349 ii++;
00350
00351 if (cn[i-1][j][2]>0) constrmat[ii][(i-1)*2+0]=-1.0;
00352 if (cn[i][j][2]>0) constrmat[ii][i*2+0]=1.0;
00353 ii++;
00354
00355 if (cn[i-1][j][5]>0) constrmat[ii][(i-1)*2+0]=-1.0;
00356 if (cn[i][j][5]>0) constrmat[ii][i*2+0]=1.0;
00357 ii++;
00358 }
00359 }
00360
00361 }
00362
00363
00364
00365
00366
00367
00368 void lplate::assemble_constr_dm ()
00369 {
00370 long i,j,k,l,*cnd,*cnm;
00371 double tl,tu;
00372 matrix lcm;
00373
00374 cnd = new long [2*ndofn];
00375 cnm = new long [nnmult];
00376 allocm (nnmult,2*ndofn,lcm);
00377
00378 fprintf (stdout,"\n pocet vsech multiplikatoru %ld",nmult);
00379 fprintf (stdout,"\n pocet vsech posunu %ld",ndof);
00380 fprintf (stdout,"\n pocet vsech posunu %ld",ndof*nproc);
00381
00382
00383 allocm (nmult,ndof*nproc,dcm);
00384
00385 for (i=1;i<nl;i++){
00386 for (j=0;j<nn;j++){
00387
00388 tl = thick[i-1][j];
00389 tu = thick[i][j];
00390
00391 for (k=0;k<ndofn;k++){
00392 cnd[k]=cn[i-1][j][k]+(i-1)*ndof;
00393 }
00394 for (k=0;k<ndofn;k++){
00395 if (cn[i][j][k]!=0)
00396 cnd[k+ndofn]=cn[i][j][k]+i*ndof;
00397 else
00398 cnd[k+ndofn]=cn[i][j][k];
00399 }
00400 for (k=0;k<nnmult;k++){
00401 cnm[k]=lgnodes[j].cn[i][k];
00402 }
00403
00404 fprintf (stdout,"\n\n");
00405 for (k=0;k<2*ndofn;k++){
00406 fprintf (stdout," %ld",cnd[k]);
00407 }
00408 fprintf (stdout,"X");
00409 for (k=0;k<nnmult;k++){
00410 fprintf (stdout," %ld",cnm[k]);
00411 }
00412
00413
00414 fillm (0.0,lcm);
00415
00416 lcm[0][0] = -1.0;
00417 lcm[1][1] = -1.0;
00418 lcm[2][2] = -1.0;
00419 lcm[1][3] = tl/2.0;
00420 lcm[0][4] = tl/(-2.0);
00421 lcm[3][5] = -1.0;
00422
00423 lcm[0][6] = 1.0;
00424 lcm[1][7] = 1.0;
00425 lcm[2][8] = 1.0;
00426 lcm[1][9] = tu/2.0;
00427 lcm[0][10] = tu/(-2.0);
00428 lcm[3][11] = 1.0;
00429
00430
00431 for (k=0;k<nnmult;k++){
00432 if (cnm[k]==0) continue;
00433 for (l=0;l<2*ndofn;l++){
00434 if (cnd[l]==0) continue;
00435 else dcm[cnm[k]-1][cnd[l]-1]=lcm[k][l];
00436 }
00437 }
00438
00439 }
00440 }
00441
00442 delete [] cnd;
00443 delete [] cnm;
00444 destrm (lcm);
00445
00446 }
00447
00448
00449
00450
00451
00452
00453
00454
00455
00456
00457
00458
00459 void lplate::orthonormalization ()
00460 {
00461 long i,j,k,ii,jj;
00462 double alpha,nor;
00463
00464
00465 for (i=0;i<nn*nnmult;i++){
00466 nor = constrmat[i][0]*constrmat[i][0] + constrmat[i][1]*constrmat[i][1];
00467 nor += constrmat[i][2]*constrmat[i][2] + constrmat[i][3]*constrmat[i][3];
00468 if (nor>zero){
00469 nor=sqrt(nor);
00470 constrmat[i][0]/=nor;
00471 constrmat[i][1]/=nor;
00472 constrmat[i][2]/=nor;
00473 constrmat[i][3]/=nor;
00474 }
00475 }
00476
00477 ii=nn*nnmult;
00478 for (i=1;i<nl-1;i++){
00479 jj=ii-nn*nnmult;
00480 for (j=0;j<nn*nnmult;j++){
00481 alpha=constrmat[ii][i*2+0]*constrmat[jj][i*2+0]+constrmat[ii][i*2+1]*constrmat[jj][i*2+1];
00482 nor=0.0;
00483 for (k=0;k<2*(i+2);k++){
00484 constrmat[ii][k]-=alpha*constrmat[jj][k];
00485 nor+=constrmat[ii][k]*constrmat[ii][k];
00486 }
00487 if (nor>zero){
00488 nor=sqrt(nor);
00489 for (k=0;k<2*(i+2);k++){
00490 constrmat[ii][k]/=nor;
00491 }
00492 }
00493 ii++; jj++;
00494 }
00495 }
00496 }
00497
00498
00499
00500
00501
00502
00503 void lplate::orthonormalization_dm ()
00504 {
00505 long i,j,k;
00506 double nor,alpha;
00507
00508 nor=0.0;
00509 for (i=0;i<nproc*ndof;i++){
00510 nor+=dcm[0][i]*dcm[0][i];
00511 }
00512 nor=sqrt(nor);
00513 if (nor>zero){
00514 for (i=0;i<nproc*ndof;i++){
00515 dcm[0][i]/=nor;
00516 }
00517 }
00518
00519 for (i=0;i<nmult;i++){
00520 for (j=i+1;j<nmult;j++){
00521 alpha=0.0;
00522 for (k=0;k<nproc*ndof;k++){
00523 alpha+=dcm[j][k]*dcm[i][k];
00524 }
00525 nor=0.0;
00526 for (k=0;k<nproc*ndof;k++){
00527 dcm[j][k]-=alpha*dcm[i][k];
00528 nor+=dcm[j][k]*dcm[j][k];
00529 }
00530 nor=sqrt(nor);
00531 if (nor>zero){
00532 for (k=0;k<nproc*ndof;k++){
00533 dcm[j][k]/=nor;
00534 }
00535 }
00536 }
00537 }
00538
00539 }
00540
00541
00542
00543
00544
00545
00546
00547
00548
00549
00550
00551 void lplate::locglob (double *gv,double *lv,long ns)
00552 {
00553 long i,j,li,ii,cnl,cn1,cn2;
00554 double cm1,cm2,lv1,lv2;
00555
00556 li=ns-1;
00557 if (li<0) li=0;
00558 ii=li*nn*nnmult;
00559 for (i=li;i<nl-1;i++){
00560 for (j=0;j<nn;j++){
00561
00562 cnl=lgnodes[j].cn[i+1][0];
00563 if (cnl==0){ ii++; }
00564 else{
00565 cm1=constrmat[ii][ns*2+0];
00566 cm2=constrmat[ii][ns*2+1];
00567 cn1=cn[ns][j][0];
00568 if (cn1>0) lv1=lv[cn1-1];
00569 else lv1=0.0;
00570 cn2=cn[ns][j][4];
00571 if (cn2>0) lv2=lv[cn2-1];
00572 else lv2=0.0;
00573 gv[cnl-1]+=cm1*lv1+cm2*lv2;
00574 ii++;
00575 }
00576
00577 cnl=lgnodes[j].cn[i+1][1];
00578 if (cnl==0){ ii++; }
00579 else{
00580 cm1=constrmat[ii][ns*2+0];
00581 cm2=constrmat[ii][ns*2+1];
00582 cn1=cn[ns][j][1];
00583 if (cn1>0) lv1=lv[cn1-1];
00584 else lv1=0.0;
00585 cn2=cn[ns][j][3];
00586 if (cn2>0) lv2=lv[cn2-1];
00587 else lv2=0.0;
00588 gv[cnl-1]+=cm1*lv1+cm2*lv2;
00589 ii++;
00590 }
00591
00592
00593 cnl=lgnodes[j].cn[i+1][2];
00594 if (cnl==0){ ii++; }
00595 else{
00596 cm1=constrmat[ii][ns*2+0];
00597 cn1=cn[ns][j][2];
00598 if (cn1>0) lv1=lv[cn1-1];
00599 else lv1=0.0;
00600 gv[cnl-1]+=cm1*lv1;
00601 ii++;
00602 }
00603
00604 cnl=lgnodes[j].cn[i+1][3];
00605 if (cnl==0){ ii++; }
00606 else{
00607 cm1=constrmat[ii][ns*2+0];
00608 cn1=cn[ns][j][5];
00609 if (cn1>0) lv1=lv[cn1-1];
00610 else lv1=0.0;
00611 gv[cnl-1]+=cm1*lv1;
00612 ii++;
00613 }
00614 }
00615 }
00616 }
00617
00618
00619
00620
00621
00622
00623
00624
00625
00626
00627
00628 void lplate::globloc (double *gv,double *lv,long ns)
00629 {
00630 long i,j,li,ii,cnl,cn1,cn2;
00631 double cm1,cm2,cgv;
00632
00633 li=ns-1;
00634 if (li<0) li=0;
00635 ii=li*nn*nnmult;
00636 for (i=li;i<nl-1;i++){
00637 for (j=0;j<nn;j++){
00638
00639 cnl=lgnodes[j].cn[i+1][0];
00640 if (cnl==0){ ii++; }
00641 else{
00642 cm1=constrmat[ii][ns*2+0];
00643 cm2=constrmat[ii][ns*2+1];
00644 cgv=gv[cnl-1];
00645
00646 cn1=cn[ns][j][0];
00647 if (cn1>0) lv[cn1-1]+=cgv*cm1;
00648 cn2=cn[ns][j][4];
00649 if (cn2>0) lv[cn2-1]+=cgv*cm2;
00650 ii++;
00651 }
00652
00653 cnl=lgnodes[j].cn[i+1][1];
00654 if (cnl==0){ ii++; }
00655 else{
00656 cm1=constrmat[ii][ns*2+0];
00657 cm2=constrmat[ii][ns*2+1];
00658 cgv=gv[cnl-1];
00659
00660 cn1=cn[ns][j][1];
00661 if (cn1>0) lv[cn1-1]+=cgv*cm1;
00662 cn2=cn[ns][j][3];
00663 if (cn2>0) lv[cn2-1]+=cgv*cm2;
00664 ii++;
00665 }
00666
00667 cnl=lgnodes[j].cn[i+1][2];
00668 if (cnl==0){ ii++; }
00669 else{
00670 cm1=constrmat[ii][ns*2+0];
00671 cgv=gv[cnl-1];
00672
00673 cn1=cn[ns][j][2];
00674 if (cn1>0) lv[cn1-1]+=cgv*cm1;
00675 ii++;
00676 }
00677
00678 cnl=lgnodes[j].cn[i+1][3];
00679 if (cnl==0){ ii++; }
00680 else{
00681 cm1=constrmat[ii][ns*2+0];
00682 cgv=gv[cnl-1];
00683
00684 cn1=cn[ns][j][5];
00685 if (cn1>0) lv[cn1-1]+=cgv*cm1;
00686 ii++;
00687 }
00688 }
00689 }
00690
00691 }
00692
00693
00694
00695
00696
00697
00698
00699
00700 void lplate::printconstrmat (FILE *out)
00701 {
00702 long i,j;
00703
00704 fprintf (out,"\n\n\n Constraint matrix \n");
00705 for (i=0;i<(nl-1)*nn*4;i++){
00706 fprintf (out,"\n row %4ld ",i);
00707 for (j=0;j<nl*2;j++){
00708 fprintf (out," %6.3e",constrmat[i][j]);
00709
00710 }
00711 }
00712 }
00713
00714
00715
00716
00717
00718
00719
00720
00721 void lplate::printconstrmat_dm (FILE *out)
00722 {
00723 long i,j;
00724
00725 fprintf (out,"\n\n\n Constraint matrix \n");
00726 for (i=0;i<nmult;i++){
00727 fprintf (out,"\n row %4ld ",i);
00728 for (j=0;j<nproc*ndof;j++){
00729 fprintf (out," %6.3e",dcm[i][j]);
00730
00731 }
00732 }
00733 }
00734
00735
00736
00737
00738
00739
00740
00741
00742
00743
00744
00745
00746
00747 void lplate::cg (gmatrix *gm,double *w,double *rhs,long *domproc,FILE *out)
00748 {
00749 long i,j;
00750 double nom,denom,alpha,beta;
00751 double *d,*dd,*r,*p,*pp;
00752 MPI_Status stat;
00753
00754 fprintf (stdout,"\n ndof %ld",ndof);
00755 fprintf (stdout,"\n maxndof %ld",maxndof);
00756 fprintf (stdout,"\n nmult %ld",nmult);
00757
00758
00759
00760 dd = new double [maxndof+1];
00761
00762 pp = new double [maxndof+1];
00763
00764 if (myrank==0){
00765 r = new double [nmult];
00766 d = new double [nmult];
00767 p = new double [nmult];
00768
00769 nullvr (w,nmult);
00770
00771
00772 for (j=1;j<nproc;j++){
00773 nullvr (dd,maxndof+1);
00774 globloc (w,dd,domproc[j]);
00775 MPI_Send (dd,maxndof,MPI_DOUBLE,j,ndom,MPI_COMM_WORLD);
00776 }
00777 nullvr (dd,maxndof+1);
00778 globloc (w,dd,domproc[0]);
00779
00780
00781 for (j=0;j<ndof;j++){
00782 dd[j]=rhs[j]-dd[j];
00783 }
00784
00785
00786 gm->back_substitution (pp,dd);
00787
00788
00789 nullvr (r,nmult);
00790 locglob (r,pp,domproc[0]);
00791
00792 for (j=1;j<nproc;j++){
00793 MPI_Recv(pp,maxndof,MPI_DOUBLE,MPI_ANY_SOURCE,MPI_ANY_TAG,MPI_COMM_WORLD,&stat);
00794 locglob (r,pp,domproc[stat.MPI_TAG]);
00795 }
00796
00797
00798
00799 for (i=0;i<nmult;i++){
00800 d[i]=r[i];
00801 }
00802
00803
00804 nom = ss (r,r,nmult);
00805
00806
00807
00808
00809
00810 for (i=0;i<nicg;i++){
00811
00812
00813
00814
00815
00816
00817 for (j=1;j<nproc;j++){
00818 nullvr (dd,maxndof+1);
00819 globloc (d,dd,domproc[j]);
00820 MPI_Send (dd,maxndof+1,MPI_DOUBLE,j,ndom,MPI_COMM_WORLD);
00821 }
00822
00823 nullvr (dd,maxndof+1);
00824 globloc (d,dd,domproc[0]);
00825
00826
00827 gm->back_substitution (pp,dd);
00828
00829
00830
00831 nullvr (p,nmult);
00832 locglob (p,pp,domproc[0]);
00833
00834 for (j=1;j<nproc;j++){
00835 MPI_Recv(pp,maxndof+1,MPI_DOUBLE,MPI_ANY_SOURCE,MPI_ANY_TAG,MPI_COMM_WORLD,&stat);
00836 locglob (p,pp,domproc[stat.MPI_TAG]);
00837 }
00838
00839
00840 denom = ss (d,p,nmult);
00841
00842 if (fabs(denom)<zero){
00843 fprintf (stderr,"\n\n zero denominator in modified conjugate gradient method (%s, line %d).\n",__FILE__,__LINE__);
00844
00845 if (i!=nicg-1){
00846 dd[maxndof]=1.0;
00847 for (j=1;j<nproc;j++){
00848 MPI_Send (dd,maxndof+1,MPI_DOUBLE,j,0,MPI_COMM_WORLD);
00849 }
00850 }
00851 break;
00852 }
00853
00854
00855
00856
00857 alpha = nom/denom;
00858
00859
00860
00861
00862 for (j=0;j<nmult;j++){
00863 w[j]+=alpha*d[j];
00864 r[j]-=alpha*p[j];
00865 }
00866
00867 denom = nom;
00868 if (fabs(denom)<zero){
00869 fprintf (stderr,"\n\n zero denominator of beta in function (%s, line %d).\n",__FILE__,__LINE__);
00870
00871 if (i!=nicg-1){
00872 dd[maxndof]=1.0;
00873 for (j=1;j<nproc;j++){
00874 MPI_Send (dd,maxndof+1,MPI_DOUBLE,j,myrank,MPI_COMM_WORLD);
00875 }
00876 }
00877 break;
00878 }
00879
00880 nom = ss (r,r,nmult);
00881
00882 fprintf (stdout,"\n iteration %4ld norm r %e",i,sqrt(nom));
00883 fprintf (out,"\n iteration %4ld norm r %e",i,sqrt(nom));
00884
00885 if (sqrt(nom)<errcg){
00886 if (i!=nicg-1){
00887 dd[maxndof]=1.0;
00888 for (j=1;j<nproc;j++){
00889 MPI_Send (dd,maxndof+1,MPI_DOUBLE,j,myrank,MPI_COMM_WORLD);
00890 }
00891 }
00892 break;
00893 }
00894
00895
00896
00897 beta = nom/denom;
00898
00899
00900 for (j=0;j<nmult;j++){
00901 d[j]=beta*d[j]+r[j];
00902 }
00903
00904 }
00905
00906
00907 anicg=i; aerrcg=nom;
00908
00909 fprintf (out,"\n\n\n\n kontrola Lagrangeovych multiplikatoru \n");
00910 for (i=0;i<nmult;i++){
00911 fprintf (out,"\n lagr. mult %4ld %e",i,w[i]);
00912 }
00913
00914 }
00915 else{
00916
00917 MPI_Recv (dd,maxndof,MPI_DOUBLE,MPI_ANY_SOURCE,MPI_ANY_TAG,MPI_COMM_WORLD,&stat);
00918
00919
00920 for (j=0;j<ndof;j++){
00921 dd[j]=rhs[j]-dd[j];
00922 }
00923
00924
00925 gm->back_substitution (pp,dd);
00926
00927 MPI_Send (pp,maxndof,MPI_DOUBLE,0,ndom,MPI_COMM_WORLD);
00928
00929
00930 for (i=0;i<nicg;i++){
00931
00932
00933 MPI_Recv (dd,maxndof+1,MPI_DOUBLE,MPI_ANY_SOURCE,MPI_ANY_TAG,MPI_COMM_WORLD,&stat);
00934 if (dd[maxndof]>0.5){ break; }
00935
00936
00937 gm->back_substitution (pp,dd);
00938
00939
00940 MPI_Send (pp,maxndof+1,MPI_DOUBLE,0,myrank,MPI_COMM_WORLD);
00941
00942 }
00943
00944 }
00945
00946 if (myrank==0){
00947 delete [] r;
00948 delete [] d;
00949 delete [] p;
00950 }
00951
00952
00953 delete [] dd;
00954 delete [] pp;
00955
00956 }
00957
00958
00959
00960
00961
00962
00963
00964
00965
00966
00967
00968
00969 void lplate::nodaldisplacements (gmatrix *gm,double *lhs,double *rhs,double *w,long *domproc)
00970 {
00971 long j;
00972 double *dd;
00973 MPI_Status stat;
00974
00975
00976 dd = new double [maxndof+1];
00977
00978 if (myrank==0){
00979
00980 for (j=1;j<nproc;j++){
00981 nullvr (dd,maxndof+1);
00982 globloc (w,dd,domproc[j]);
00983 MPI_Send (dd,maxndof,MPI_DOUBLE,j,ndom,MPI_COMM_WORLD);
00984 }
00985
00986 nullvr (dd,maxndof+1);
00987 globloc (w,dd,domproc[0]);
00988
00989
00990 for (j=0;j<ndof;j++){
00991 dd[j]=rhs[j]-dd[j];
00992 }
00993
00994
00995 gm->back_substitution (lhs,dd);
00996
00997 }
00998 else{
00999
01000 MPI_Recv (dd,maxndof,MPI_DOUBLE,MPI_ANY_SOURCE,MPI_ANY_TAG,MPI_COMM_WORLD,&stat);
01001
01002
01003 for (j=0;j<ndof;j++){
01004 dd[j]=rhs[j]-dd[j];
01005 }
01006
01007
01008 gm->back_substitution (lhs,dd);
01009
01010 }
01011
01012 delete [] dd;
01013 }
01014
01015
01016
01017
01018
01019
01020
01021
01022
01023
01024 void lplate::solve_system (gtopology *top,gmatrix *gm,
01025 long *domproc,double *lhs,double *rhs,FILE *out)
01026 {
01027 double *w;
01028
01029 fprintf (out,"\n\n kontrola prave strany\n");
01030 long i;
01031 for (i=0;i<ndof;i++){
01032 fprintf (out," %e",rhs[i]);
01033 }
01034
01035
01036
01037 if (myrank==0)
01038 w = new double [nmult];
01039
01040
01041 gm->decompose_matrix ();
01042
01043
01044
01045
01046
01047 if (myrank==0) orthonormalization ();
01048
01049
01050
01051
01052
01053
01054 cg (gm,w,rhs,domproc,out);
01055
01056
01057 nodaldisplacements (gm,lhs,rhs,w,domproc);
01058
01059 if (myrank==0)
01060 delete [] w;
01061
01062 }