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