00001 #include "partop.h"
00002 #include <string.h>
00003 #include <math.h>
00004 #include <time.h>
00005 #include "../GEFEL/galias.h"
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015 partop::partop(int np,int mr,long nd,meshdescription meshd)
00016 {
00017
00018 nproc=np;
00019
00020 myrank=mr;
00021
00022 ndom=nd;
00023
00024
00025 md = meshd;
00026
00027
00028 nn=0;
00029
00030 ne=0;
00031
00032 tnnp=0;
00033
00034 tnbn=0;
00035
00036 nn=0;
00037
00038 maxnn=0;
00039
00040 nbn=0;
00041
00042 maxnbn = 0;
00043
00044 maxndof=0;
00045
00046
00047
00048 nodmultip = NULL;
00049
00050 lnbn = NULL;
00051
00052 lnin = NULL;
00053
00054 lgnbn = NULL;
00055
00056 dofmultip = NULL;
00057
00058 nodeidentif=NULL;
00059
00060 if (myrank==0){
00061 nnsd = NULL;
00062 nbnd = NULL;
00063 multip = NULL;
00064 allnodes = NULL;
00065 cnbn = NULL;
00066 nbdofnd = NULL;
00067 nbdofd = NULL;
00068 lbcn = NULL;
00069 bnmultip = NULL;
00070 llnbn = NULL;
00071 ldn = NULL;
00072
00073
00074 nalldof=NULL;
00075
00076 gnn = NULL;
00077 pgcn = NULL;
00078 gcnbn = NULL;
00079 }
00080
00081 }
00082
00083
00084
00085
00086
00087
00088 partop::~partop()
00089 {
00090 long i;
00091
00092 delete [] nodmultip;
00093 delete [] lnbn;
00094
00095 if (myrank==0){
00096 delete [] nnsd;
00097 delete [] nbnd;
00098 delete [] multip;
00099 delete [] nbdofd;
00100 delete [] bnmultip;
00101
00102
00103 for (i=0;i<nproc;i++){
00104 delete [] allnodes[i];
00105 }
00106 delete [] allnodes;
00107
00108 for (i=0;i<nproc;i++){
00109 delete [] cnbn[i];
00110 }
00111 delete [] cnbn;
00112
00113 for (i=0;i<nproc;i++){
00114 delete [] nbdofnd[i];
00115 }
00116 delete [] nbdofnd;
00117
00118 for (i=0;i<tnbn;i++){
00119 delete [] llnbn[i];
00120 }
00121 delete [] llnbn;
00122
00123 for (i=0;i<tnbn;i++){
00124 delete [] ldn[i];
00125 }
00126 delete [] ldn;
00127
00128
00129 }
00130 }
00131
00132
00133
00134
00135
00136
00137
00138
00139
00140
00141 void partop::initiation (gtopology *top,long *ltg)
00142 {
00143 long i;
00144
00145 nn = top->nn;
00146 ne = top->ne;
00147
00148
00149 switch (md){
00150 case all_nodes:{
00151 for (i=0;i<nn;i++){
00152 top->gnodes[i].ai = ltg[i];
00153 }
00154 break;
00155 }
00156 case bound_nodes:{
00157 for (i=0;i<nn;i++){
00158 top->gnodes[i].ai = ltg[i];
00159 }
00160 break;
00161 }
00162 case neg_bound_nodes:{
00163 for (i=0;i<nn;i++){
00164 top->gnodes[i].ai = abs(ltg[i]);
00165 }
00166 break;
00167 }
00168 default:{
00169 fprintf (stderr,"\n\n unknown type of mesh description is required");
00170 fprintf (stderr,"\n in function initiation (file %s, line %d)\n",__FILE__,__LINE__);
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 void partop::numbers_of_all_nodes_on_subdomains (long *domproc,FILE *out)
00196 {
00197 long i,j,k;
00198 MPI_Status stat;
00199
00200
00201 maxnn=nn;
00202
00203
00204
00205
00206 if (myrank==0){
00207
00208 if (nnsd!=NULL)
00209 delete [] nnsd;
00210 nnsd = new long [nproc];
00211
00212
00213 j=domproc[0];
00214 nnsd[j]=maxnn;
00215
00216 for (i=1;i<nproc;i++){
00217 MPI_Recv (&k,1,MPI_LONG,MPI_ANY_SOURCE,MPI_ANY_TAG,MPI_COMM_WORLD,&stat);
00218
00219
00220 j=domproc[stat.MPI_TAG];
00221 if (maxnn<k) maxnn=k;
00222 nnsd[j]=k;
00223 }
00224
00225 for (i=1;i<nproc;i++){
00226 MPI_Send (&maxnn,1,MPI_LONG,i,myrank,MPI_COMM_WORLD);
00227 }
00228 }
00229
00230
00231
00232 else{
00233 MPI_Send (&nn,1,MPI_LONG,0,myrank,MPI_COMM_WORLD);
00234 MPI_Recv (&maxnn,1,MPI_LONG,MPI_ANY_SOURCE,MPI_ANY_TAG,MPI_COMM_WORLD,&stat);
00235 }
00236 MPI_Barrier (MPI_COMM_WORLD);
00237
00238
00239 if (myrank==0){
00240 fprintf (out,"\n\n\n numbers of nodes on each subdomain, maxnn %ld\n\n",maxnn);
00241 for (i=0;i<nproc;i++){
00242 fprintf (out,"\n %ld",nnsd[i]);
00243 }
00244 }
00245
00246 }
00247
00248
00249
00250
00251
00252
00253
00254
00255
00256
00257
00258
00259
00260
00261
00262
00263 void partop::numbers_of_all_dofs_on_subdomains (gtopology *top,long *domproc,FILE *out)
00264 {
00265 long i,j,k;
00266 MPI_Status stat;
00267
00268
00269 maxndof=0;
00270 for (i=0;i<nn;i++){
00271 maxndof+=top->give_ndofn (i);
00272 }
00273
00274
00275
00276
00277 if (myrank==0){
00278
00279 if (nalldof!=NULL)
00280 delete [] nalldof;
00281 nalldof = new long [nproc];
00282
00283
00284 j=domproc[0];
00285 nalldof[j]=maxndof;
00286
00287 for (i=1;i<nproc;i++){
00288 MPI_Recv (&k,1,MPI_LONG,MPI_ANY_SOURCE,MPI_ANY_TAG,MPI_COMM_WORLD,&stat);
00289
00290
00291 j=domproc[stat.MPI_TAG];
00292 if (maxndof<k) maxndof=k;
00293 nalldof[j]=k;
00294 }
00295
00296 for (i=1;i<nproc;i++){
00297 MPI_Send (&maxndof,1,MPI_LONG,i,myrank,MPI_COMM_WORLD);
00298 }
00299 }
00300
00301
00302
00303 else{
00304 MPI_Send (&maxndof,1,MPI_LONG,0,myrank,MPI_COMM_WORLD);
00305 MPI_Recv (&maxndof,1,MPI_LONG,MPI_ANY_SOURCE,MPI_ANY_TAG,MPI_COMM_WORLD,&stat);
00306 }
00307 MPI_Barrier (MPI_COMM_WORLD);
00308
00309
00310 if (myrank==0){
00311 fprintf (out,"\n\n\n numbers of all DOFs on each subdomain, maxndof %ld\n\n",maxndof);
00312 for (i=0;i<nproc;i++){
00313 fprintf (out,"\n %ld",nalldof[i]);
00314 }
00315 }
00316
00317 }
00318
00319
00320
00321
00322
00323
00324
00325
00326
00327
00328
00329
00330
00331
00332
00333
00334
00335
00336
00337
00338
00339
00340
00341
00342
00343
00344 void partop::compute_multiplicity (long *ltg,long *domproc,FILE *out)
00345 {
00346 long i,j,k,a,buffsize;
00347 long *buff;
00348 MPI_Status stat;
00349
00350 if (maxnn==0){
00351 fprintf (stderr,"\n\n maximum number of nodes on one subdomain is equal to zero,\n");
00352 fprintf (stderr,"\n function numbers_of_all_nodes_on_subdomains has to be called first,\n");
00353 abort ();
00354 }
00355
00356 switch (md){
00357 case bound_nodes:{
00358
00359 buffsize=2;
00360 buff = new long [buffsize];
00361
00362
00363 maxnbn=0;
00364 tnbn=0;
00365 for (i=0;i<nn;i++){
00366 if (tnbn<ltg[i])
00367 tnbn=ltg[i];
00368 if (ltg[i]>-1)
00369 maxnbn++;;
00370 }
00371
00372 buff[0]=tnbn;
00373 buff[1]=maxnbn;
00374
00375 if (myrank==0){
00376 for (i=1;i<nproc;i++){
00377 MPI_Recv (buff,buffsize,MPI_LONG,MPI_ANY_SOURCE,MPI_ANY_TAG,MPI_COMM_WORLD,&stat);
00378 if (tnbn<buff[0])
00379 tnbn=buff[0];
00380 if (maxnbn<buff[1])
00381 maxnbn=buff[1];
00382 }
00383 tnbn++;
00384
00385 buff[0]=tnbn;
00386 buff[1]=maxnbn;
00387
00388 for (i=1;i<nproc;i++){
00389 MPI_Send (buff,buffsize,MPI_LONG,i,myrank,MPI_COMM_WORLD);
00390 }
00391 }
00392 else{
00393 MPI_Send (buff,buffsize,MPI_LONG,0,myrank,MPI_COMM_WORLD);
00394 MPI_Recv (buff,buffsize,MPI_LONG,MPI_ANY_SOURCE,MPI_ANY_TAG,MPI_COMM_WORLD,&stat);
00395 }
00396 MPI_Barrier (MPI_COMM_WORLD);
00397
00398 tnbn=buff[0];
00399 maxnbn=buff[1];
00400
00401 delete [] buff;
00402
00403 buffsize = maxnbn;
00404 buff = new long [buffsize];
00405 for (i=0;i<maxnbn;i++){
00406 buff[i]=-1;
00407 }
00408
00409 j=0;
00410 for(i=0;i<nn;i++){
00411 if (ltg[i]>-1){
00412 buff[j] = ltg[i];
00413 j++;
00414 }
00415 }
00416
00417
00418 if(myrank == 0){
00419 if (multip!=NULL){
00420 delete [] multip;
00421 }
00422 multip = new long [tnbn];
00423 for(i=0;i<tnbn;i++){
00424 multip[i]=0;
00425 }
00426
00427
00428 for(j=0;j<maxnbn;j++){
00429 if (buff[j]>-1){
00430 multip[buff[j]]++;
00431 }
00432 }
00433
00434
00435 for (i=1;i<nproc;i++){
00436 MPI_Recv (buff,buffsize,MPI_LONG,MPI_ANY_SOURCE,MPI_ANY_TAG,MPI_COMM_WORLD,&stat);
00437 for(j=0;j<maxnbn;j++){
00438 if (buff[j]>-1){
00439 multip[buff[j]]++;
00440 }
00441 }
00442 }
00443
00444 }
00445 else{
00446 MPI_Send (buff,buffsize,MPI_LONG,0,myrank,MPI_COMM_WORLD);
00447 }
00448 MPI_Barrier (MPI_COMM_WORLD);
00449 delete [] buff;
00450
00451 buffsize=tnbn;
00452 buff = new long [buffsize];
00453
00454 if (myrank==0){
00455 for (i=0;i<tnbn;i++){
00456 buff[i]=multip[i];
00457 }
00458
00459 for (i=1;i<nproc;i++){
00460 MPI_Send (buff,buffsize,MPI_LONG,i,myrank,MPI_COMM_WORLD);
00461 }
00462 }
00463 else{
00464 MPI_Recv (buff,buffsize,MPI_LONG,MPI_ANY_SOURCE,MPI_ANY_TAG,MPI_COMM_WORLD,&stat);
00465 }
00466 MPI_Barrier (MPI_COMM_WORLD);
00467
00468 if (nodmultip!=NULL)
00469 delete [] nodmultip;
00470 nodmultip = new long [nn];
00471
00472 for (i=0;i<nn;i++){
00473 if (ltg[i]>-1){
00474 nodmultip[i]=buff[ltg[i]];
00475 }
00476 else{
00477 nodmultip[i]=1;
00478 }
00479 }
00480
00481 delete [] buff;
00482
00483 break;
00484 }
00485
00486
00487
00488
00489 case all_nodes:
00490 case neg_bound_nodes:{
00491
00492 buffsize=maxnn+1;
00493 buff = new long [buffsize];
00494 tnnp=0;
00495 for (i=0;i<nn;i++){
00496 buff[i]=ltg[i];
00497 if (tnnp<ltg[i])
00498 tnnp=ltg[i];
00499 }
00500 buff[maxnn]=tnnp;
00501
00502
00503
00504
00505 if (myrank==0){
00506 if (allnodes != NULL){
00507 for (i=0;i<nproc;i++){
00508 delete [] allnodes[i];
00509 }
00510 delete [] allnodes;
00511 }
00512 allnodes = new long* [nproc];
00513 for (i=0;i<nproc;i++){
00514 allnodes[i] = new long [nnsd[i]];
00515 }
00516
00517
00518 k=domproc[0];
00519 for (j=0;j<nnsd[k];j++){
00520 allnodes[k][j]=buff[j];
00521 }
00522 if (tnnp<buff[maxnn])
00523 tnnp=buff[maxnn];
00524
00525
00526 for (i=1;i<nproc;i++){
00527 MPI_Recv (buff,buffsize,MPI_LONG,MPI_ANY_SOURCE,MPI_ANY_TAG,MPI_COMM_WORLD,&stat);
00528
00529 k=domproc[stat.MPI_TAG];
00530 for (j=0;j<nnsd[k];j++){
00531 allnodes[k][j]=buff[j];
00532 }
00533 if (tnnp<buff[maxnn])
00534 tnnp=buff[maxnn];
00535 }
00536 }
00537
00538
00539
00540
00541 else{
00542 MPI_Send (buff,buffsize,MPI_LONG,0,myrank,MPI_COMM_WORLD);
00543 }
00544
00545 if (myrank==0){
00546
00547
00548
00549
00550 tnnp++;
00551
00552 fprintf (out,"\n\n\n total number of nodes on whole problem %ld\n",tnnp);
00553 if (multip != NULL)
00554 delete [] multip;
00555 multip = new long [tnnp];
00556 for (i=0;i<tnnp;i++){
00557 multip[i]=0;
00558 }
00559
00560
00561 for (i=0;i<nproc;i++){
00562 for (j=0;j<nnsd[i];j++){
00563 multip[allnodes[i][j]]++;
00564 }
00565 }
00566 }
00567
00568
00569 if (myrank==0){
00570 for (i=1;i<nproc;i++){
00571 k=domproc[i];
00572 for (j=0;j<nnsd[k];j++){
00573 buff[j]=multip[allnodes[k][j]];
00574 }
00575 MPI_Send (buff,buffsize,MPI_LONG,i,myrank,MPI_COMM_WORLD);
00576 }
00577
00578 k=domproc[0];
00579 for (j=0;j<nnsd[k];j++){
00580 buff[j]=multip[allnodes[k][j]];
00581 }
00582 }
00583 else{
00584 MPI_Recv (buff,buffsize,MPI_LONG,MPI_ANY_SOURCE,MPI_ANY_TAG,MPI_COMM_WORLD,&stat);
00585 }
00586
00587
00588 if (nodmultip != NULL)
00589 delete [] nodmultip;
00590 nodmultip = new long [nn];
00591 for (i=0;i<nn;i++){
00592 nodmultip[i] = buff[i];
00593 }
00594
00595 delete [] buff;
00596 break;
00597 }
00598 default:{
00599 fprintf (stderr,"\n\n unknown type of mesh description is required in");
00600 fprintf (stderr,"\n function compute_multiplicity (file %s, line %d)\n",__FILE__,__LINE__);
00601 }
00602 }
00603 MPI_Barrier (MPI_COMM_WORLD);
00604
00605
00606
00607
00608
00609
00610
00611
00612
00613
00614
00615
00616
00617
00618
00619
00620
00621
00622
00623
00624 }
00625
00626
00627
00628
00629
00630
00631
00632
00633 void partop::dof_multiplicity (gtopology *top,FILE *out)
00634 {
00635 long i,j,k,nm;
00636 long ndofn;
00637
00638 ndof = 0;
00639 for (i=0;i<nn;i++){
00640 ndofn=top->give_ndofn (i);
00641 for (j=0;j<ndofn;j++){
00642 k=top->give_dof (i,j);
00643 if (k>0){
00644 ndof++;
00645 }
00646 }
00647 }
00648
00649
00650
00651
00652 dofmultip = new long [ndof];
00653
00654 for (i=0;i<nn;i++){
00655 ndofn=top->give_ndofn (i);
00656 nm=nodmultip[i];
00657 for (j=0;j<ndofn;j++){
00658 k=top->give_dof (i,j);
00659 if (k>0){
00660 dofmultip[k-1]=nm;
00661 }
00662 }
00663 }
00664
00665
00666
00667
00668
00669
00670
00671 }
00672
00673
00674
00675
00676
00677
00678
00679
00680
00681
00682
00683
00684
00685
00686
00687
00688
00689
00690
00691
00692
00693
00694
00695
00696
00697
00698
00699
00700
00701 void partop::find_boundary_nodes (long *ltg,long *domproc,FILE *out)
00702 {
00703 long i,j,k,l,buffsize;
00704 long *buff;
00705 MPI_Status stat;
00706
00707
00708
00709
00710
00711
00712 switch (md){
00713 case bound_nodes:{
00714
00715 nbn=0; maxnbn=0;
00716 for (i=0;i<nn;i++){
00717 if (ltg[i]>-1) nbn++;
00718 if (maxnbn<ltg[i]) maxnbn=ltg[i];
00719 }
00720
00721 buffsize=2;
00722 buff = new long [buffsize];
00723 buff[0]=nbn;
00724 buff[1]=maxnbn;
00725
00726
00727
00728
00729
00730 if (myrank==0){
00731
00732 if (nbnd != NULL)
00733 delete [] nbnd;
00734 nbnd = new long [nproc];
00735
00736
00737 j=domproc[0];
00738 if (maxnbn<buff[0]) maxnbn=buff[0];
00739 if (tnbn<buff[1]) tnbn=buff[1];
00740 nbnd[j]=buff[0];
00741
00742 for (i=1;i<nproc;i++){
00743 MPI_Recv (buff,buffsize,MPI_LONG,MPI_ANY_SOURCE,MPI_ANY_TAG,MPI_COMM_WORLD,&stat);
00744
00745
00746 j=domproc[stat.MPI_TAG];
00747 if (maxnbn<buff[0]) maxnbn=buff[0];
00748 if (tnbn<buff[1]) tnbn=buff[1];
00749 nbnd[j]=buff[0];
00750
00751 }
00752
00753
00754 tnbn++;
00755
00756 for (i=1;i<nproc;i++){
00757 MPI_Send (&maxnbn,1,MPI_LONG,i,myrank,MPI_COMM_WORLD);
00758 }
00759 }
00760
00761
00762
00763 else{
00764 MPI_Send (buff,buffsize,MPI_LONG,0,myrank,MPI_COMM_WORLD);
00765 MPI_Recv (&maxnbn,1,MPI_LONG,MPI_ANY_SOURCE,MPI_ANY_TAG,MPI_COMM_WORLD,&stat);
00766 }
00767
00768 delete [] buff;
00769 fprintf (out,"\n\n maxnbn is %ld",maxnbn);
00770 break;
00771 }
00772
00773 case all_nodes:
00774 case neg_bound_nodes:{
00775
00776 if (myrank==0){
00777 if (multip==NULL){
00778 fprintf (stderr,"\n\n array multip is not allocated in function find_boundary_nodes,\n");
00779 fprintf (stderr,"\n function compute_multiplicity has to be called first (file %s, line %d),\n",__FILE__,__LINE__);
00780 }
00781 }
00782
00783 if (myrank==0){
00784
00785 tnbn=0;
00786 for (i=0;i<tnnp;i++){
00787 if (multip[i]!=1) tnbn++;
00788 }
00789
00790
00791 if (nbnd != NULL)
00792 delete [] nbnd;
00793 nbnd = new long [nproc];
00794
00795 maxnbn=0;
00796 for (i=0;i<nproc;i++){
00797 k=0;
00798 for (j=0;j<nnsd[i];j++){
00799 if (multip[allnodes[i][j]]>1) k++;
00800 }
00801 if (maxnbn<k) maxnbn=k;
00802 nbnd[i]=k;
00803 }
00804
00805 for (i=1;i<nproc;i++){
00806 MPI_Send (&maxnbn,1,MPI_LONG,i,myrank,MPI_COMM_WORLD);
00807 }
00808 }
00809 else{
00810 MPI_Recv (&maxnbn,1,MPI_LONG,MPI_ANY_SOURCE,MPI_ANY_TAG,MPI_COMM_WORLD,&stat);
00811 }
00812
00813
00814
00815 nbn=0;
00816 for (i=0;i<nn;i++){
00817 if (nodmultip[i]>1)
00818 nbn++;
00819 }
00820
00821 break;
00822 }
00823 default:{
00824 fprintf (stderr,"\n\n unknown type of mesh description is required in");
00825 fprintf (stderr,"\n function find_boundary_nodes (file %s, line %d)\n",__FILE__,__LINE__);
00826 }
00827 }
00828
00829
00830
00831 buffsize=maxnbn;
00832 buff = new long [buffsize];
00833
00834 switch (md){
00835 case bound_nodes:{
00836
00837
00838 if (lgnbn!=NULL)
00839 delete [] lgnbn;
00840 lgnbn = new long [nbn];
00841
00842 j=0;
00843 for (i=0;i<nn;i++){
00844 if (ltg[i]>-1){
00845 lgnbn[j]=ltg[i];
00846 j++;
00847 }
00848 }
00849 break;
00850 }
00851 case all_nodes:
00852 case neg_bound_nodes:{
00853
00854 if (myrank==0){
00855 for (i=1;i<nproc;i++){
00856 k=domproc[i];
00857 l=0;
00858 for (j=0;j<nnsd[k];j++){
00859 if (multip[allnodes[k][j]]>1){
00860 buff[l]=allnodes[k][j];
00861 l++;
00862 }
00863 }
00864
00865 MPI_Send (buff,buffsize,MPI_LONG,i,myrank,MPI_COMM_WORLD);
00866 }
00867
00868 k=domproc[0];
00869 l=0;
00870 for (j=0;j<nnsd[k];j++){
00871 if (multip[allnodes[k][j]]>1){
00872 buff[l]=allnodes[k][j];
00873 l++;
00874 }
00875 }
00876
00877 }
00878 else{
00879 MPI_Recv (buff,buffsize,MPI_LONG,MPI_ANY_SOURCE,MPI_ANY_TAG,MPI_COMM_WORLD,&stat);
00880 }
00881
00882
00883 if (lgnbn!=NULL)
00884 delete [] lgnbn;
00885 lgnbn = new long [nbn];
00886
00887 for (i=0;i<nbn;i++){
00888 lgnbn[i]=buff[i];
00889 }
00890
00891 break;
00892 }
00893 default:{
00894 fprintf (stderr,"\n\n unknown type of mesh description is required in");
00895 fprintf (stderr,"\n function find_boundary_nodes (file %s, line %d)\n",__FILE__,__LINE__);
00896 }
00897 }
00898
00899
00900 delete [] buff;
00901
00902
00903
00904
00905
00906
00907 if (lnbn != NULL)
00908 delete [] lnbn;
00909 lnbn = new long [nbn];
00910
00911 switch (md){
00912 case all_nodes:
00913 case neg_bound_nodes:{
00914 j=0;
00915 for (i=0;i<nn;i++){
00916 if (nodmultip[i]>1){
00917 lnbn[j]=i;
00918 j++;
00919 }
00920 }
00921
00922 break;
00923 }
00924
00925 case bound_nodes:{
00926 j=0;
00927 for (i=0;i<nn;i++){
00928 if (ltg[i]>-1){
00929 lnbn[j]=i;
00930 j++;
00931 }
00932 }
00933 break;
00934 }
00935 default:{
00936 fprintf (stderr,"\n\n unknown type of mesh description is required in");
00937 fprintf (stderr,"\n function find_boundary_nodes (file %s, line %d)\n",__FILE__,__LINE__);
00938 }
00939 }
00940
00941
00942
00943 nin = nn-nbn;
00944
00945
00946 if (lnin != NULL)
00947 delete [] lnin;
00948 lnin = new long [nin];
00949
00950 buff = new long [nn];
00951 for (i=0;i<nn;i++){
00952 buff[i]=0;
00953 }
00954
00955 for (i=0;i<nbn;i++){
00956 j=lnbn[i];
00957 buff[j]=1;
00958 }
00959
00960 j=0;
00961 for (i=0;i<nn;i++){
00962 if (buff[i]==0){
00963 lnin[j]=i;
00964 j++;
00965 }
00966 }
00967
00968 delete [] buff;
00969
00970 MPI_Barrier (MPI_COMM_WORLD);
00971
00972
00973 if (myrank==0){
00974 fprintf (out,"\n\n\n numbers of boundary nodes on each subdomain \n\n");
00975 for (i=0;i<nproc;i++){
00976 fprintf (out,"\n %ld",nbnd[i]);
00977 }
00978 }
00979
00980 fprintf (out,"\n\n\n numbers of boundary nodes on subdomain \n\n");
00981 fprintf (out,"\n %ld",nbn);
00982 fprintf (out,"\n\n\n local numbers of boundary nodes on subdomain");
00983 for (i=0;i<nbn;i++){
00984 fprintf (out,"\n %ld %ld",i,lnbn[i]+1);
00985 }
00986
00987 fprintf (out,"\n\n\n numbers of internal nodes on subdomain \n\n");
00988 fprintf (out,"\n %ld",nin);
00989 fprintf (out,"\n\n\n local numbers of internal nodes on subdomain");
00990 for (i=0;i<nin;i++){
00991 fprintf (out,"\n %ld %ld",i,lnin[i]+1);
00992 }
00993
00994
00995 }
00996
00997
00998
00999
01000
01001
01002 void partop::rewrite_ltg (long *ltg)
01003 {
01004 long i;
01005
01006 for (i=0;i<nn;i++){
01007 ltg[i]=-1;
01008 }
01009 for (i=0;i<nbn;i++){
01010 ltg[lnbn[i]]=lgnbn[i];
01011 }
01012 }
01013
01014
01015
01016
01017
01018
01019
01020
01021
01022
01023
01024
01025
01026
01027
01028
01029
01030
01031
01032
01033
01034 void partop::boundary_nodes_on_master (gtopology *top,long *domproc,FILE *out)
01035 {
01036 long i,j,k,m,buffsize;
01037 long *buff;
01038 MPI_Status stat;
01039
01040 if (myrank==0){
01041 if (nbnd==NULL){
01042 fprintf (stderr,"\n\n array containing numbers of boundary nodes on subdomains is not allocated,\n");
01043 fprintf (stderr,"\n function find_boundary_nodes has to be called first (file %s, line %d),\n",__FILE__,__LINE__);
01044 }
01045 }
01046
01047
01048 if (myrank==0){
01049 if (cnbn != NULL){
01050 for (i=0;i<nproc;i++){
01051 delete [] cnbn[i];
01052 }
01053 delete [] cnbn;
01054 }
01055 cnbn = new long* [nproc];
01056 for (i=0;i<nproc;i++){
01057 cnbn[i] = new long [nbnd[i]];
01058 }
01059 }
01060
01061 switch (md){
01062 case all_nodes:
01063 case neg_bound_nodes:{
01064
01065 if (myrank==0){
01066 if (gcnbn != NULL)
01067 delete [] gcnbn;
01068 gcnbn = new long [tnnp];
01069
01070 j=0;
01071 for (i=0;i<tnnp;i++){
01072 if (multip[i]>1){
01073 gcnbn[i]=j;
01074 j++;
01075 }
01076 else{
01077 gcnbn[i]=-1;
01078 }
01079 }
01080
01081
01082
01083
01084
01085
01086
01087
01088
01089 if (j!=tnbn){
01090 fprintf (stderr,"\n\n total number of boundary nodes computed in function boundary_nodes_on_master");
01091 fprintf (stderr,"\n differs from the number computed in function find_boundary_nodes (file %s, line %d),\n",__FILE__,__LINE__);
01092 }
01093
01094 for (i=0;i<nproc;i++){
01095 m=0;
01096 for (j=0;j<nnsd[i];j++){
01097 k=allnodes[i][j];
01098 if (gcnbn[k]>-1){
01099 cnbn[i][m]=gcnbn[k];
01100 m++;
01101 }
01102 }
01103 }
01104 delete [] gcnbn;
01105 gcnbn = NULL;
01106 }
01107
01108 break;
01109 }
01110
01111 case bound_nodes:{
01112
01113 if (maxnbn==0){
01114 fprintf (stderr,"\n\n maximum number of boundary nodes on one subdomain is equal to zero,");
01115 fprintf (stderr,"\n function find_boundary_nodes should be called first (file %s, line %d),\n",__FILE__,__LINE__);
01116 }
01117
01118 buffsize=maxnbn;
01119 buff = new long [buffsize];
01120
01121
01122 j=0;
01123 for (i=0;i<nn;i++){
01124 if (top->gnodes[i].ai>-1){
01125 buff[j]=top->gnodes[i].ai;
01126 j++;
01127 }
01128 }
01129
01130 if (myrank==0){
01131
01132
01133 k=domproc[0];
01134 for (j=0;j<nbnd[k];j++){
01135 cnbn[k][j]=buff[j];
01136 }
01137
01138 for (i=1;i<nproc;i++){
01139 MPI_Recv (buff,buffsize,MPI_LONG,MPI_ANY_SOURCE,MPI_ANY_TAG,MPI_COMM_WORLD,&stat);
01140
01141
01142 k=domproc[stat.MPI_TAG];
01143 for (j=0;j<nbnd[k];j++){
01144 cnbn[k][j]=buff[j];
01145 }
01146
01147 }
01148
01149 }
01150 else{
01151 MPI_Send (buff,buffsize,MPI_LONG,0,myrank,MPI_COMM_WORLD);
01152 }
01153
01154 delete [] buff;
01155 break;
01156 }
01157 default:{
01158 fprintf (stderr,"\n\n unknown type of mesh description is required in");
01159 fprintf (stderr,"\n function boundary_nodes_on_master (file %s, line %d)\n",__FILE__,__LINE__);
01160 }
01161 }
01162
01163
01164 MPI_Barrier (MPI_COMM_WORLD);
01165
01166
01167 if (myrank==0){
01168 fprintf (out,"\n\n\n\n coarse numbering of boundary nodes - cnbn\n\n");
01169 for (i=0;i<nproc;i++){
01170 fprintf (out,"\n\n domain number %ld\n",i+1);
01171 for (j=0;j<nbnd[i];j++){
01172 fprintf (out,"%ld %ld\n",j,cnbn[i][j]+1);
01173 }
01174 }
01175 }
01176
01177
01178 }
01179
01180
01181
01182
01183
01184
01185
01186
01187
01188
01189
01190
01191
01192
01193 void partop::coarse_local_nodes ()
01194 {
01195 long i,j,k,m;
01196
01197 if (myrank==0){
01198
01199 bnmultip = new long [tnbn];
01200 for (i=0;i<tnbn;i++){
01201 bnmultip[i]=0;
01202 }
01203
01204 if (md == all_nodes){
01205 j=0;
01206 for (i=0;i<tnnp;i++){
01207 if (multip[i]>1){
01208 bnmultip[j]=multip[i];
01209 j++;
01210 }
01211 }
01212 }
01213
01214 if (md == bound_nodes){
01215 for (i=0;i<nproc;i++){
01216 for (j=0;j<nbnd[i];j++){
01217 bnmultip[cnbn[i][j]]++;
01218 }
01219 }
01220 }
01221
01222 llnbn = new long* [tnbn];
01223 ldn = new long* [tnbn];
01224 for (i=0;i<tnbn;i++){
01225 llnbn[i] = new long [bnmultip[i]];
01226 ldn[i] = new long [bnmultip[i]];
01227 bnmultip[i]=0;
01228 }
01229
01230 for (i=0;i<nproc;i++){
01231 for (j=0;j<nbnd[i];j++){
01232 k=cnbn[i][j];
01233 m=bnmultip[k];
01234 llnbn[k][m]=j;
01235 ldn[k][m]=i;
01236 bnmultip[k]++;
01237 }
01238 }
01239 }
01240 }
01241
01242
01243
01244
01245
01246
01247
01248
01249 void partop::sort_nodes (FILE *out)
01250 {
01251 long i,j,k,m,n,min;
01252
01253 if (myrank==0){
01254 for (i=0;i<tnbn;i++){
01255 for (j=0;j<bnmultip[i];j++){
01256 min=nproc;
01257 for (k=j;k<bnmultip[i];k++){
01258 if (ldn[i][k]<min){
01259 min=ldn[i][k];
01260 m=k;
01261 }
01262 }
01263 n=ldn[i][j];
01264 ldn[i][j]=ldn[i][m];
01265 ldn[i][m]=n;
01266
01267 n=llnbn[i][j];
01268 llnbn[i][j]=llnbn[i][m];
01269 llnbn[i][m]=n;
01270 }
01271 }
01272 }
01273
01274 if (myrank==0){
01275 if(md == all_nodes){
01276 fprintf (out,"\n\n Multiplicity of all nodes on whole problem\n");
01277 for (i=0;i<tnnp;i++){
01278 fprintf (out,"%ld %ld\n",i+1,multip[i]);
01279 }
01280 }
01281 fprintf (out,"\n\n Multiplicity of all boundary nodes on whole problem\n");
01282 for (i=0;i<tnbn;i++){
01283 fprintf (out,"%ld %ld\n",i+1,bnmultip[i]);
01284 }
01285
01286
01287
01288
01289
01290
01291
01292
01293
01294
01295
01296 }
01297
01298 }
01299
01300
01301
01302
01303
01304
01305
01306
01307
01308
01309
01310
01311
01312
01313
01314
01315
01316
01317
01318
01319
01320
01321
01322 void partop::number_of_bdofs_on_nodes (gtopology *top,long *domproc,FILE *out)
01323 {
01324 long i,j,k,buffsize;
01325 long *buff;
01326 MPI_Status stat;
01327
01328 if (lnbn==NULL){
01329 fprintf (stderr,"\n\n array containing local numbers of boundary nodes is not allocated,\n");
01330 fprintf (stderr,"\n function find_boundary_nodes has to be called first (file %s, line %d),\n",__FILE__,__LINE__);
01331 }
01332
01333 buffsize = maxnbn;
01334 buff = new long [buffsize];
01335
01336
01337 for (i=0;i<nbn;i++){
01338 j=lnbn[i];
01339 buff[i]=top->give_ndofn (j);
01340 }
01341
01342 if (myrank==0){
01343 if (nbdofnd != NULL){
01344 for (i=0;i<nproc;i++){
01345 delete [] nbdofnd[i];
01346 }
01347 delete [] nbdofnd;
01348 }
01349 nbdofnd = new long* [nproc];
01350 for (i=0;i<nproc;i++){
01351 nbdofnd[i] = new long [nbnd[i]];
01352 }
01353
01354
01355 k=domproc[0];
01356 for (j=0;j<nbnd[k];j++){
01357 nbdofnd[k][j]=buff[j];
01358 }
01359
01360 for (i=1;i<nproc;i++){
01361 MPI_Recv (buff,buffsize,MPI_LONG,MPI_ANY_SOURCE,MPI_ANY_TAG,MPI_COMM_WORLD,&stat);
01362
01363
01364 k=domproc[stat.MPI_TAG];
01365 for (j=0;j<nbnd[k];j++){
01366 nbdofnd[k][j]=buff[j];
01367 }
01368
01369 }
01370 }
01371 else{
01372 MPI_Send (buff,buffsize,MPI_LONG,0,myrank,MPI_COMM_WORLD);
01373 }
01374 MPI_Barrier (MPI_COMM_WORLD);
01375
01376 delete [] buff;
01377
01378
01379 if (myrank==0){
01380 fprintf (out,"\n\n\n numbers of DOFs at boundary nodes");
01381 for (i=0;i<nproc;i++){
01382 fprintf (out,"\n\n domain %ld",i);
01383 for (j=0;j<nbnd[i];j++){
01384 fprintf (out,"\n %ld %ld",j+1,nbdofnd[i][j]);
01385 }
01386 }
01387 }
01388
01389 }
01390
01391
01392
01393
01394
01395
01396
01397
01398
01399
01400
01401
01402
01403
01404
01405
01406
01407
01408
01409
01410
01411
01412
01413
01414
01415
01416 void partop::code_numbers_on_master (gtopology *top,long *domproc,FILE *out)
01417 {
01418 long i,j,k,m,l,ndofn,buffsize;
01419 long *buff;
01420 MPI_Status stat;
01421
01422 if (myrank==0){
01423 if (nbdofnd==NULL){
01424 fprintf (stderr,"\n\n array ndofnd is not allocated in function code_numbers_on_master,\n");
01425 fprintf (stderr,"\n function number_of_dofs_on_nodes has to be called first (file %s, line %d),\n",__FILE__,__LINE__);
01426 }
01427 }
01428
01429
01430 if (myrank==0){
01431 if (nbdofd != NULL)
01432 delete [] nbdofd;
01433 nbdofd = new long [nproc];
01434
01435 maxnbdof=0;
01436 for (i=0;i<nproc;i++){
01437 nbdofd[i]=0;
01438 for (j=0;j<nbnd[i];j++){
01439 nbdofd[i]+=nbdofnd[i][j];
01440 }
01441 if (maxnbdof<nbdofd[i])
01442 maxnbdof=nbdofd[i];
01443 }
01444
01445 for (i=1;i<nproc;i++){
01446 MPI_Send (&maxnbdof,1,MPI_LONG,i,myrank,MPI_COMM_WORLD);
01447 }
01448 }
01449 else{
01450 MPI_Recv (&maxnbdof,1,MPI_LONG,MPI_ANY_SOURCE,MPI_ANY_TAG,MPI_COMM_WORLD,&stat);
01451 }
01452
01453 buffsize = maxnbdof;
01454 buff = new long [buffsize];
01455
01456 m=0;
01457 for (i=0;i<nbn;i++){
01458 k=lnbn[i];
01459 ndofn = top->give_ndofn (k);
01460 for (j=0;j<ndofn;j++){
01461 buff[m]=top->give_dof (k,j);
01462 m++;
01463 }
01464 }
01465
01466 if (myrank==0){
01467 if (lbcn != NULL){
01468 for (i=0;i<nproc;i++){
01469 for (j=0;j<nbnd[i];j++){
01470 delete [] lbcn[i][j];
01471 }
01472 delete [] lbcn[i];
01473 }
01474 delete [] lbcn;
01475 }
01476 lbcn = new long** [nproc];
01477 for (i=0;i<nproc;i++){
01478 lbcn[i] = new long* [nbnd[i]];
01479 for (j=0;j<nbnd[i];j++){
01480 lbcn[i][j] = new long [nbdofnd[i][j]];
01481 }
01482 }
01483
01484
01485 l=0;
01486 m=domproc[0];
01487 for (j=0;j<nbnd[m];j++){
01488 for (k=0;k<nbdofnd[m][j];k++){
01489 lbcn[m][j][k]=buff[l];
01490 l++;
01491 }
01492 }
01493
01494 for (i=1;i<nproc;i++){
01495 MPI_Recv (buff,buffsize,MPI_LONG,MPI_ANY_SOURCE,MPI_ANY_TAG,MPI_COMM_WORLD,&stat);
01496
01497
01498 l=0;
01499 m=domproc[stat.MPI_TAG];
01500 for (j=0;j<nbnd[m];j++){
01501 for (k=0;k<nbdofnd[m][j];k++){
01502 lbcn[m][j][k]=buff[l];
01503 l++;
01504 }
01505 }
01506
01507 }
01508 }
01509 else{
01510 MPI_Send (buff,buffsize,MPI_LONG,0,myrank,MPI_COMM_WORLD);
01511 }
01512
01513
01514 if (myrank==0){
01515 fprintf (out,"\n\n\n kontrola pole lbcn \n\n");
01516 for (i=0;i<nproc;i++){
01517 fprintf (out,"\n\n domena %ld",i);
01518 for (j=0;j<nbnd[i];j++){
01519 fprintf (out,"\n boundary uzel %ld",j);
01520 for (k=0;k<nbdofnd[i][j];k++){
01521 fprintf (out," %ld",lbcn[i][j][k]);
01522 }
01523 fprintf (out,"\n");
01524 }
01525 }
01526 }
01527
01528
01529
01530 MPI_Barrier (MPI_COMM_WORLD);
01531
01532 }
01533
01534
01535
01536
01537
01538
01539
01540
01541 long partop::give_whole_elemtype(gtopology *top)
01542 {
01543 long i,t;
01544 gelemtype get;
01545 long *elem;
01546 long type = 0;
01547 elem = new long[ne];
01548
01549 for(i = 0; i < ne; i++){
01550 get = top->give_elem_type (i);
01551 switch(get){
01552 case linbar:{ t= 1; break; }
01553 case quadbar:{ t= 2; break; }
01554 case lintriag:{ t=21; break; }
01555 case quadtriag:{ t=22; break; }
01556 case linquad:{ t=25; break; }
01557 case quadquad:{ t=26; break; }
01558 case lintetra:{ t=41; break; }
01559 case quadtetra:{ t=42; break; }
01560 case linhexa:{ t=45; break; }
01561 case quadhexa:{ t=46; break; }
01562 default:{
01563 fprintf (stderr,"\n\n unknown element type is required in function give_whole_elemtype (file %s, line %d).\n",__FILE__,__LINE__);
01564 }
01565 }
01566 elem[i] = t;
01567 }
01568
01569 for(i = 1; i < ne; i++){
01570 if(elem[i] != elem[i-1]){
01571 fprintf (stderr,"\n\n inhomogeneous mesh, function identification_node can not be used\n");
01572 }
01573 else{
01574 type = elem[i];
01575 }
01576 }
01577 return type;
01578 delete []elem;
01579
01580 }
01581
01582
01583
01584
01585
01586
01587
01588
01589 long partop::give_whole_dim(gtopology *top)
01590 {
01591
01592 gelemtype get;
01593 long *elem;
01594 long dim,i,t;
01595 dim = -1;
01596
01597 elem = new long[ne];
01598
01599 for(i = 0; i < ne; i++){
01600 get = top->give_elem_type (i);
01601 switch(get){
01602 case linbar:{ t = 1; break; }
01603 case quadbar:{ t = 1; break; }
01604 case lintriag:{ t = 2; break; }
01605 case quadtriag:{ t = 2; break; }
01606 case linquad:{ t = 2; break; }
01607 case quadquad:{ t = 2; break; }
01608 case lintetra:{ t = 3; break; }
01609 case quadtetra:{ t = 3; break; }
01610 case linhexa:{ t = 3; break; }
01611 case quadhexa:{ t = 3; break; }
01612 default:{
01613 fprintf (stderr,"\n\n unknown element type is required in function give_whole_dim (file %s, line %d).\n",__FILE__,__LINE__);
01614 }
01615 }
01616 elem[i] = t;
01617 }
01618
01619
01620 for(i = 1; i < ne; i++){
01621 if(elem[i] != elem[i-1]){
01622 fprintf (stderr,"\n\n different spatial dimension is found in the mesh, function identification_node can not be used\n");
01623 }
01624 else{
01625 dim = elem[i];
01626 }
01627 }
01628
01629 return(dim);
01630
01631 delete []elem;
01632
01633 }
01634
01635
01636 void partop::identification_node(gtopology *top,long *ltg,long *domproc, FILE *out)
01637 {
01638 long dim;
01639 dim = give_whole_dim(top);
01640
01641 fprintf(out,"\n\n\nDimension of problem is %ld\n",dim);
01642 switch(dim){
01643 case -1:
01644 fprintf (stderr,"\n\n Spatial dimension is not detected (file %s, line %d).\n",__FILE__,__LINE__);
01645 break;
01646 case 1:
01647 break;
01648 case 2: identification_node_2d (top,ltg,domproc,out);
01649 break;
01650 case 3: identification_node_3d (top,ltg,domproc,out);
01651 break;
01652 default:
01653 fprintf (stderr,"\n unknown dimension is used in function identification_node (file %s, line %d)\n",__FILE__,__LINE__);
01654 break;
01655 }
01656 }
01657
01658
01659
01660
01661
01662
01663
01664
01665
01666
01667
01668
01669
01670
01671
01672
01673 void partop::identification_node_2d (gtopology *top, long *ltg,long *domproc, FILE *out)
01674 {
01675 long i,j,k,l,m,n,a,b;
01676
01677 if (nodeidentif != NULL)
01678 delete []nodeidentif;
01679 nodeidentif = new long [nn];
01680
01681 for (i=0;i<nn;i++){
01682 nodeidentif[i]=1;
01683 }
01684
01685
01686
01687 top->adjacnodes_edge (out);
01688
01689
01690
01691
01692
01693 for (i = 0; i < nn; i++){
01694 if (nodmultip[i] == 1)
01695
01696 nodeidentif[i] = 1;
01697 if(nodmultip[i] > 2){
01698
01699 nodeidentif[i] = 3;
01700 }
01701
01702
01703
01704
01705
01706
01707
01708
01709
01710
01711
01712
01713
01714
01715
01716
01717
01718
01719
01720
01721 if(nodmultip[i] == 2){
01722 m = 0;
01723 b = 0;
01724
01725 for (j = 0;j < top->nadjacnodesedge[i]; j++){
01726 l = 0;
01727 for(k = 0; k < top->n_nned[i][j]; k++){
01728 a = top->adjacnodesedge[i][j][k];
01729 if(a != i){
01730
01731 if(nodmultip[a] == 2){
01732 l++;
01733 }
01734 if(nodmultip[a] > 2){
01735 b++;
01736 }
01737 }
01738 }
01739
01740 if(l == 1){
01741 m++;
01742 }
01743 }
01744
01745 if(m == 1 && b == 0){
01746
01747 nodeidentif[i] = 3;
01748 }
01749 else{
01750
01751 nodeidentif[i] = 2;
01752 }
01753 }
01754 if(nodmultip[i] > 2){
01755 m = 0;
01756 b = 0;
01757 n = 0;
01758
01759 for (j = 0;j < top->nadjacnodesedge[i]; j++){
01760 l = 0;
01761 for(k = 0; k < top->n_nned[i][j]; k++){
01762 a = top->adjacnodesedge[i][j][k];
01763 if(a != i){
01764
01765 if(nodmultip[i] == nodmultip[a]){
01766 l++;
01767 }
01768 if( nodmultip[i] < nodmultip[a] ){
01769 b++;
01770 }
01771 if(nodmultip[i] > nodmultip[a]){
01772 n++;
01773 }
01774 }
01775 }
01776
01777 if( l == 1){
01778 m++;
01779 }
01780 }
01781
01782
01783 if((m == 1 && b == 0) || n == top -> nadjacnodesedge[i] ){
01784
01785 nodeidentif[i] = 3;
01786 }
01787 else{
01788
01789 nodeidentif[i] = 2;
01790 }
01791
01792
01793 }
01794 }
01795
01796
01797
01798 for (i = 0; i < nn; i++){
01799 if(nodeidentif[i] == 0) fprintf(out,"\n node number %ld is not identification",i+1);
01800 if(nodeidentif[i] == 1) fprintf(out,"\n node number %ld is internal node",i+1);
01801 if(nodeidentif[i] == 2) fprintf(out,"\n node number %ld is edge node",i+1);
01802 if(nodeidentif[i] == 3) fprintf(out,"\n node number %ld is vertex node",i+1);
01803 }
01804
01805
01806 for(i = 0; i < nn; i++){
01807 for(j = 0; j < top->nadjacnodesedge[i]; j++){
01808 delete []top->adjacnodesedge[i][j];
01809 }
01810 delete []top->adjacnodesedge[i];
01811 delete []top->n_nned[i];
01812 }
01813 delete []top->adjacnodesedge;
01814 delete []top->n_nned;
01815 delete []top->nadjacnodesedge;
01816
01817 }
01818
01819
01820
01821
01822
01823
01824
01825
01826
01827
01828
01829
01830
01831
01832
01833
01834
01835 void partop::identification_node_3d(gtopology *top,long *ltg,long *domproc, FILE *out)
01836 {
01837 long type;
01838
01839 type = give_whole_elemtype(top);
01840 switch(type){
01841 case 45:
01842 fprintf(out,"\n\n\nUsed element into problem is hexahedron\n\n",type);
01843 identification_node_hexa (top,ltg,domproc,out);
01844 break;
01845 case 46:
01846 fprintf(out,"\n\n\nUsed element into problem is hexahedron\n\n",type);
01847 identification_node_hexa (top,ltg,domproc,out);
01848 break;
01849 case 41:
01850 fprintf(out,"\n\n\nUsed element into problem is tetrahedron\n\n",type);
01851 identification_node_tetra (top, domproc,out);
01852 break;
01853 case 42:
01854 fprintf(out,"\n\n\nUsed element into problem is tetrahedron\n\n",type);
01855 identification_node_tetra (top, domproc,out);
01856 break;
01857 default:
01858 fprintf (stderr,"\n\n unknown element type is required in function identification_node_3d (file %s, line %d).\n",__FILE__,__LINE__);
01859 break;
01860 }
01861 }
01862
01863
01864
01865
01866
01867
01868
01869
01870
01871
01872
01873
01874
01875
01876
01877
01878
01879
01880
01881
01882 void partop::identification_node_hexa (gtopology *top,long *ltg,long *domproc, FILE *out)
01883 {
01884 long a,i,j,k,l,m,b,n;
01885
01886 if (nodeidentif != NULL)
01887 delete []nodeidentif;
01888 nodeidentif = new long [nn];
01889
01890 for (i = 0; i < nn; i++)
01891 nodeidentif[i] = 0;
01892
01893
01894 top->adjacnodes_edge (out);
01895
01896
01897
01898
01899
01900
01901 for (i=0;i<nn;i++){
01902
01903 if(nodmultip[i] == 1){
01904
01905 nodeidentif[i] = 1;
01906 }
01907
01908
01909
01910
01911
01912
01913
01914
01915
01916
01917
01918
01919
01920
01921
01922
01923
01924
01925
01926
01927
01928
01929
01930
01931
01932
01933
01934
01935
01936
01937
01938
01939
01940
01941
01942
01943
01944
01945
01946
01947
01948
01949
01950
01951
01952
01953
01954
01955
01956
01957
01958
01959
01960
01961
01962 if(nodmultip[i] == 2){
01963 m = 0;
01964 b = 0;
01965
01966 for (j = 0;j < top->nadjacnodesedge[i]; j++){
01967 l = 0;
01968 for(k = 0; k < top->n_nned[i][j]; k++){
01969 a = top->adjacnodesedge[i][j][k];
01970 if(a != i){
01971
01972 if(nodmultip[a] == 2){
01973 l++;
01974 }
01975 if(nodmultip[a] > 2){
01976 b++;
01977 }
01978 }
01979 }
01980
01981 if(l == 1){
01982 m++;
01983 }
01984 }
01985
01986 if((m == 2 && top->nadjacnodesedge[i] == 3) && b == 0){
01987
01988 nodeidentif[i] = 3;
01989 }
01990 else{
01991
01992 nodeidentif[i] = 2;
01993 }
01994 }
01995 if(nodmultip[i] >= 3){
01996 m = 0;
01997 b = 0;
01998 n = 0;
01999
02000 for (j = 0;j < top->nadjacnodesedge[i]; j++){
02001 l = 0;
02002 for(k = 0; k < top->n_nned[i][j]; k++){
02003 a = top->adjacnodesedge[i][j][k];
02004 if(a != i){
02005
02006 if(nodmultip[i] == nodmultip[a]){
02007 l++;
02008 }
02009 if( nodmultip[i] < nodmultip[a] ){
02010 b++;
02011 }
02012 if(nodmultip[i] > nodmultip[a]){
02013 n++;
02014 }
02015 }
02016 }
02017 if( l == 1){
02018 m++;
02019 }
02020 }
02021
02022 if((m == 1 && b == 0) || n == top -> nadjacnodesedge[i] ){
02023
02024 nodeidentif[i] = 3;
02025 }
02026 else{
02027
02028 nodeidentif[i] = 2;
02029 }
02030 }
02031
02032 }
02033
02034 for (i=0;i<nn;i++){
02035 if(nodeidentif[i]==0) fprintf(out,"\n node number %ld is not identification",i+1);
02036 if(nodeidentif[i]==1) fprintf(out,"\n node number %ld is internal node",i+1);
02037 if(nodeidentif[i]==2) fprintf(out,"\n node number %ld is edge node",i+1);
02038 if(nodeidentif[i]==4) fprintf(out,"\n node number %ld is face node",i+1);
02039 if(nodeidentif[i]==3) fprintf(out,"\n node number %ld is vertex node",i+1);
02040 }
02041
02042
02043 for(i = 0; i < nn; i++){
02044 for(j = 0; j < top->nadjacnodesedge[i]; j++){
02045 delete []top->adjacnodesedge[i][j];
02046 }
02047 delete []top->adjacnodesedge[i];
02048 delete []top->n_nned[i];
02049 }
02050 delete []top->adjacnodesedge;
02051 delete []top->n_nned;
02052 delete []top->nadjacnodesedge;
02053
02054
02055 }
02056
02057
02058
02059
02060
02061
02062
02063
02064
02065
02066
02067
02068
02069
02070
02071
02072
02073
02074
02075
02076 void partop::identification_node_tetra (gtopology *top, long *domproc, FILE *out)
02077 {
02078 long i,j,k,m,l,a,b,n;
02079
02080 if (nodeidentif != NULL)
02081 delete [] nodeidentif;
02082 nodeidentif = new long [nn];
02083
02084 for (i=0;i<nn;i++) nodeidentif[i] = 0;
02085
02086 top->adjacnodes_edge(out);
02087
02088
02089
02090
02091
02092 for (i=0;i<nn;i++){
02093 if(nodmultip[i] == 1){
02094
02095 nodeidentif[i]=1;
02096 }
02097
02098
02099 if(nodmultip[i] == 2){
02100 m = 0;
02101 b = 0;
02102 n = 0;
02103
02104 for (j = 0;j < top->nadjacnodesedge[i]; j++){
02105 l = 0;
02106 for(k = 0; k < top->n_nned[i][j]; k++){
02107 a = top->adjacnodesedge[i][j][k];
02108 if(a != i){
02109
02110 if(nodmultip[a] == 2){
02111 l++;
02112 }
02113 if(nodmultip[a] > 2){
02114 b++;
02115 }
02116 if(nodmultip[a] < 2){
02117 n++;
02118 }
02119 }
02120 }
02121
02122 if(l == 1){
02123 m++;
02124 }
02125 }
02126
02127 if(m == 2 && b == 0 ){
02128
02129 nodeidentif[i] = 3;
02130 }
02131 else{
02132
02133 nodeidentif[i] = 2;
02134 }
02135 }
02136 if(nodmultip[i] >= 3){
02137 m = 0;
02138 b = 0;
02139 n = 0;
02140
02141 for (j = 0;j < top->nadjacnodesedge[i]; j++){
02142 l = 0;
02143 for(k = 0; k < top->n_nned[i][j]; k++){
02144 a = top->adjacnodesedge[i][j][k];
02145 if(a != i){
02146
02147 if(nodmultip[i] == nodmultip[a]){
02148 l++;
02149 }
02150 if( nodmultip[i] < nodmultip[a] ){
02151 b++;
02152 }
02153 if(nodmultip[i] > nodmultip[a]){
02154 n++;
02155 }
02156 }
02157 }
02158
02159 if( l == 1){
02160 m++;
02161 }
02162 }
02163
02164
02165 if((m == 1 && b == 0) || n == top -> nadjacnodesedge[i] ){
02166
02167 nodeidentif[i] = 3;
02168 }
02169 else{
02170
02171 nodeidentif[i] = 2;
02172 }
02173 }
02174 }
02175
02176
02177 for (i=0;i<nn;i++){
02178 if(nodeidentif[i]==0) fprintf(out,"\n node number %ld is not identification",i+1);
02179 if(nodeidentif[i]==1) fprintf(out,"\n node number %ld is internal node",i+1);
02180 if(nodeidentif[i]==2) fprintf(out,"\n node number %ld is edge node",i+1);
02181 if(nodeidentif[i]==4) fprintf(out,"\n node number %ld is face node",i+1);
02182 if(nodeidentif[i]==3) fprintf(out,"\n node number %ld is vertex node",i+1);
02183 }
02184
02185
02186 for(i = 0; i < nn; i++){
02187 for(j = 0; j < top->nadjacnodesedge[i]; j++){
02188 delete []top->adjacnodesedge[i][j];
02189 }
02190 delete []top->adjacnodesedge[i];
02191 delete []top->n_nned[i];
02192 }
02193 delete []top->adjacnodesedge;
02194 delete []top->n_nned;
02195 delete []top->nadjacnodesedge;
02196
02197 }
02198
02199 void partop::control_numbers_of_vertices (gtopology *top,FILE *out)
02200 {
02201 long i,j,k,l,m,a,b,dim;
02202 long n,nom,denom,nc;
02203 long *vertex,**vcomb,*line;
02204 double u1,u2,u3,v1,v2,v3,x1,x2,x3,y1,y2,y3,z1,z2,z3,scal,u,v,cosalpha;
02205
02206
02207
02208
02209 n = 0;
02210 for (i=0;i<nn;i++){
02211 if(nodeidentif[i]==3) n++;
02212 }
02213
02214 fprintf(out,"\n\n\nNumber of vertices is %ld\n\n",n);
02215
02216 dim = give_whole_dim(top);
02217
02218 if(dim == 3){
02219 if( n < 3){
02220 fprintf (stderr,"Number of vertices is < 3\n");
02221 fprintf (stderr,"Problem in function control_numbers_of_vertices %s, line %d).\n",__FILE__,__LINE__);
02222 }
02223 }
02224 if(dim == 2){
02225 if( n < 2){
02226 fprintf (stderr,"Number of vertices is < 3\n");
02227 fprintf (stderr,"Problem in function control_numbers_of_vertices %s, line %d).\n",__FILE__,__LINE__);
02228 }
02229 }
02230
02231
02232 vertex = new long[n];
02233 m = 0;
02234 for (i=0;i<nn;i++){
02235 if(nodeidentif[i]==3){
02236 vertex[m] = i;
02237
02238 m++;
02239 }
02240 }
02241
02242
02243
02244
02245
02246
02247
02248 if (n > 3){
02249
02250
02251
02252
02253 nom = 1;
02254 for(i = 1; i < n+1; i++){
02255 nom*=i;
02256 }
02257 denom = 1;
02258 for(i = 1; i < n-2; i++){
02259 denom*=i;
02260 }
02261 nc = nom/(denom*6);
02262 fprintf(out,"Pocet kombinaci k vyzkouseni je %ld\n",nc);
02263
02264
02265 vcomb = new long*[nc];
02266 for(i = 0; i < nc; i++ ){
02267 vcomb[i] = new long[3];
02268 }
02269 line = new long[nc];
02270
02271
02272 b = 0;
02273 a = n-2;
02274 for(i = 0; i < n-2; i++){
02275
02276
02277 x1 = top->gnodes[vertex[i]].x;
02278 y1 = top->gnodes[vertex[i]].y;
02279 z1 = top->gnodes[vertex[i]].z;
02280 l = i+1;
02281 for(j = 0; j < a; j++){
02282
02283
02284 x2 = top->gnodes[vertex[l]].x;
02285 y2 = top->gnodes[vertex[l]].y;
02286 z2 = top->gnodes[vertex[l]].z;
02287
02288 u1 = x1 - x2;
02289
02290 u2 = y1 - y2;
02291
02292 u3 = z1 - z2;
02293
02294 m = l+1;
02295 for(k = 0; k < a-j; k++){
02296
02297 vcomb[b][0]= vertex[i];
02298 vcomb[b][1]= vertex[l];
02299 vcomb[b][2]= vertex[m];
02300
02301 x3 = top->gnodes[vertex[m]].x;
02302 y3 = top->gnodes[vertex[m]].y;
02303 z3 = top->gnodes[vertex[m]].z;
02304
02305 v1 = x1 - x3;
02306
02307 v2 = y1 - y3;
02308
02309 v3 = z1 - z3;
02310
02311
02312 scal = u1*v1+u2*v2+u3*v3;
02313
02314
02315 u = u1*u1+u2*u2+u3*u3;
02316
02317 u = sqrt(u);
02318
02319
02320 v = v1*v1+v2*v2+v3*v3;
02321
02322 v = sqrt(v);
02323
02324
02325 cosalpha = fabs(scal/(u*v));
02326
02327 if(cosalpha == 1.0000000){
02328
02329 line[b] = 1;
02330 b++;
02331 }
02332 else{
02333 line[b] = 0;
02334 b++;
02335 }
02336 m++;
02337 }
02338 l++;
02339 }
02340 a--;
02341 }
02342
02343 b = 0;
02344 for(i = 0; i < nc; i++){
02345 if(line[i] == 0){
02346 fprintf(out,"body %ld %ld %ld nelezi na primce\n",vcomb[i][0]+1,vcomb[i][1]+1,vcomb[i][2]+1);
02347 }
02348 else{
02349 b++;
02350 fprintf(out,"body %ld %ld %ld lezi na primce\n",vcomb[i][0]+1,vcomb[i][1]+1,vcomb[i][2]+1);
02351 }
02352 }
02353 fprintf(out," %ld kombinaci lezi na primce\n",b);
02354 }
02355
02356 delete []vertex;
02357 }
02358
02359 void partop::control_vertices_on_master(gtopology *top, long *domproc, FILE *out)
02360 {
02361 long i,j,k,m,n;
02362 long buffsize;
02363 long *buff,*aux,*auxnn;
02364 long **globnodeidentif;
02365 MPI_Status stat;
02366
02367
02368
02369
02370 buffsize = maxnbn;
02371 buff = new long[buffsize];
02372
02373 j = 0;
02374 for(i = 0; i < nn; i++ ){
02375 if( nodeidentif[i] > 1){
02376 buff[j] = nodeidentif[i];
02377 j++;
02378 }
02379 }
02380
02381 if(myrank == 0){
02382
02383 globnodeidentif = new long*[nproc];
02384
02385
02386 j=domproc[0];
02387
02388
02389 globnodeidentif[j] = new long[nbnd[j]];
02390 for(i = 0; i < nbnd[0]; i++){
02391 globnodeidentif[j][i] = buff[i];
02392 }
02393
02394
02395
02396 for (k=1;k<nproc;k++){
02397 MPI_Recv (buff,buffsize,MPI_LONG,MPI_ANY_SOURCE,MPI_ANY_TAG,MPI_COMM_WORLD,&stat);
02398 j=domproc[stat.MPI_TAG];
02399 globnodeidentif[j] = new long[nbnd[j]];
02400 for(i = 0; i < nbnd[j]; i++){
02401 globnodeidentif[j][i] = buff[i];
02402 }
02403 }
02404
02405
02406
02407
02408
02409
02410
02411
02412
02413
02414 aux = new long[nproc];
02415 auxnn = new long[nproc];
02416 for(i = 0; i < tnbn; i++){
02417
02418 for(j = 0; j < nproc; j++){
02419 aux[j] = 0;
02420 auxnn[j] = -1;
02421 for(k = 0; k < nbnd[j]; k++){
02422 if(cnbn[j][k] == i){
02423 aux[j] = globnodeidentif[j][k];
02424 auxnn[j] = k;
02425 break;
02426 }
02427 }
02428 }
02429
02430
02431
02432
02433
02434
02435
02436 m = 0;
02437
02438 n = 0;
02439 for(j = 0; j < nproc; j++){
02440 if(aux[j] == 2){
02441 n++;
02442 }
02443 if(aux[j] == 3){
02444 m++;
02445 }
02446 }
02447 if( m >= n){
02448
02449 for(j = 0; j < nproc; j++){
02450 if(aux[j] != 0 && auxnn[j] != -1){
02451 globnodeidentif[j][auxnn[j]] = 3;
02452 }
02453 }
02454 }
02455 else{
02456 for(j = 0; j < nproc; j++){
02457 if(aux[j] != 0 && auxnn[j] != -1){
02458 globnodeidentif[j][auxnn[j]] = 2;
02459 }
02460 }
02461 }
02462 }
02463 delete []aux;
02464 delete []auxnn;
02465
02466
02467
02468
02469
02470
02471
02472
02473
02474
02475 for (i = 1;i < nproc; i++){
02476 for(k = 0; k < nbnd[i]; k++){
02477 buff[k] = globnodeidentif[i][k];
02478 }
02479 MPI_Send (buff,buffsize,MPI_LONG,i,myrank,MPI_COMM_WORLD);
02480 }
02481
02482
02483 for(k = 0; k < nbnd[0]; k++){
02484 buff[k] = globnodeidentif[0][k];
02485 }
02486
02487
02488 for (k=0;k<nproc;k++){
02489 delete []globnodeidentif[k];
02490 }
02491 delete []globnodeidentif;
02492
02493 }
02494 else{
02495 MPI_Send (buff,buffsize,MPI_LONG,0,myrank,MPI_COMM_WORLD);
02496 MPI_Recv (buff,buffsize,MPI_LONG,MPI_ANY_SOURCE,MPI_ANY_TAG,MPI_COMM_WORLD,&stat);
02497 }
02498 MPI_Barrier (MPI_COMM_WORLD);
02499
02500 for(i = 0; i < nbn; i++){
02501 if(buff[i] == nodeidentif[lnbn[i]]){
02502 nodeidentif[lnbn[i]] = buff[i];
02503 }
02504 else{
02505 fprintf(out,"pozor zmena\n");
02506 nodeidentif[lnbn[i]] = buff[i];
02507 }
02508 }
02509 delete []buff;
02510
02511 }
02512
02513
02514
02515
02516
02517
02518 void partop::corner_detection(gtopology *top, long *domproc, long *ltg, FILE *out)
02519 {
02520 long i,j,k;
02521 long buffsize,c,b;
02522 long *buff,*master_corner;
02523 MPI_Status stat;
02524
02525
02526 identification_node(top,ltg,domproc,out);
02527 control_numbers_of_vertices (top,out);
02528 control_vertices_on_master(top,domproc,out);
02529
02530
02531 buffsize = maxnbn;
02532 buff = new long[buffsize];
02533
02534
02535 k = 0;
02536 for(i = 0; i < nn; i++){
02537 if(nodeidentif[i] > 1){
02538 buff[k] = nodeidentif[i];
02539 k++;
02540
02541 }
02542 }
02543
02544
02545
02546 if(myrank==0){
02547
02548 master_corner = new long[tnbn];
02549
02550 for(i = 0; i < tnbn; i++ ) master_corner[i] = -1;
02551
02552
02553 j=domproc[0];
02554 for(k = 0; k < nbnd[j]; k++){
02555 master_corner[cnbn[j][k]]=buff[k];
02556 }
02557
02558 for (i=1;i<nproc;i++){
02559 MPI_Recv (buff,buffsize,MPI_LONG,MPI_ANY_SOURCE,MPI_ANY_TAG,MPI_COMM_WORLD,&stat);
02560
02561
02562 j=domproc[stat.MPI_TAG];
02563 for(k = 0; k < nbnd[j]; k++){
02564 if(master_corner[cnbn[j][k]] == -1){
02565 master_corner[cnbn[j][k]]= buff[k];
02566 }
02567 else{
02568 if(master_corner[cnbn[j][k]] != buff[k]){
02569 fprintf (out,"\n %ld %ld %ld different nodeidentif on the %ld-th domain in the %ld-th node",cnbn[j][k],buff[k],master_corner[cnbn[j][k]],j+1,k);
02570 fprintf (stderr,"different nodeidentif on the %ld-th domain in the %ld-th node\n",j,k);
02571 fprintf (stderr,"\n in function corner_detection (file %s, line %d)\n",__FILE__,__LINE__);
02572 }
02573 }
02574 }
02575 }
02576
02577
02578
02579
02580
02581
02582
02583 c = 1;
02584 b = 1;
02585 for(i = 0; i < tnbn; i++){
02586 if(master_corner[i] == 3){
02587 master_corner[i] = c*(-1);
02588 c++;
02589 }
02590 if(master_corner[i] == 2){
02591 master_corner[i] = b;
02592 b++;
02593 }
02594
02595 }
02596
02597 fprintf(out,"\n\n\n Number of corner nodes in whole problem is %ld\n",c-1);
02598 fprintf(out,"\n\n\n Number of boundary nodes in whole problem is %ld\n",b-1);
02599
02600 if((c-1)+(b-1) != tnbn){
02601 fprintf (stderr,"Problem in function corner_detection - the number of corner and edge nodes is not equal to total number of boundary nodes (file %s, line %d).\n",__FILE__,__LINE__);
02602
02603 }
02604
02605
02606 for (i = 1;i < nproc; i++){
02607 for(k = 0; k < nbnd[i]; k++){
02608 buff[k] = master_corner[cnbn[i][k]];
02609 }
02610 MPI_Send (buff,buffsize,MPI_LONG,i,myrank,MPI_COMM_WORLD);
02611 }
02612
02613
02614 for(k = 0; k < nbnd[0]; k++){
02615 buff[k] = master_corner[cnbn[0][k]];
02616 }
02617 delete []master_corner;
02618 }
02619
02620
02621 else{
02622 MPI_Send (buff,buffsize,MPI_LONG,0,myrank,MPI_COMM_WORLD);
02623 MPI_Recv (buff,buffsize,MPI_LONG,MPI_ANY_SOURCE,MPI_ANY_TAG,MPI_COMM_WORLD,&stat);
02624 }
02625 MPI_Barrier (MPI_COMM_WORLD);
02626
02627
02628
02629
02630
02631
02632
02633
02634 fprintf(out,"\n\n\n array corneridentif\n");
02635
02636 k = 0;
02637 for(i = 0; i < nn; i++){
02638 if(nodeidentif[i] == 1 || nodeidentif[i] == 0){
02639 ltg[i] = 0;
02640 }
02641 else{
02642 ltg[i] = buff[k];
02643 k++;
02644 }
02645 fprintf(out,"%ld %ld\n",i+1,ltg[i]);
02646
02647 ltg[i]--;
02648 }
02649
02650 delete []buff;
02651 delete []nodeidentif;
02652
02653 }
02654
02655
02656
02657
02658
02659
02660
02661
02662
02663
02664
02665
02666
02667 void partop::assemble_gnn(gtopology *top,long *domproc,FILE *out)
02668 {
02669 long i,j,k;
02670 long a,b,buffsize;
02671 long *buff;
02672 MPI_Status stat;
02673
02674
02675
02676
02677
02678
02679
02680
02681
02682
02683
02684
02685
02686
02687
02688
02689
02690
02691
02692
02693
02694
02695
02696
02697
02698
02699
02700
02701
02702
02703
02704
02705
02706
02707
02708
02709
02710
02711
02712
02713
02714
02715
02716
02717
02718
02719
02720
02721
02722
02723
02724
02725
02726
02727
02728
02729
02730
02731
02732
02733
02734
02735
02736
02737
02738
02739
02740
02741
02742
02743
02744
02745
02746
02747
02748
02749
02750
02751 buffsize = maxnn;
02752 buff = new long [buffsize];
02753
02754
02755 for (i=0;i<nn;i++){
02756 buff[i]=top -> give_ndofn(i);
02757 }
02758
02759 if (myrank==0){
02760
02761
02762 if (gnn == NULL)
02763 gnn = new long[tnnp+1];
02764
02765 for(i = 0; i < tnnp+1; i++)
02766 gnn[i] = -10;
02767
02768
02769 k = domproc[0];
02770 for (i = 0;i < nnsd[k]; i++){
02771 gnn[allnodes[k][i]] = buff[i];
02772 }
02773
02774 for (i= 1 ;i < nproc; i++){
02775 MPI_Recv (buff,buffsize,MPI_LONG,MPI_ANY_SOURCE,MPI_ANY_TAG,MPI_COMM_WORLD,&stat);
02776
02777
02778 k=domproc[stat.MPI_TAG];
02779 for (j = 0; j < nnsd[k]; j++){
02780 if(gnn[allnodes[k][j]] == -10){
02781 gnn[allnodes[k][j]] = buff[j];
02782 }
02783 else{
02784 if (gnn[allnodes[k][j]] != buff[j]){
02785 fprintf (stderr,"\n different numbers of DOFs on the %ld-th domain in the %ld-th node",k,j);
02786 fprintf (stderr,"\n in function assemble_gnn (file %s, line %d)\n",__FILE__,__LINE__);
02787 }
02788 }
02789 }
02790 }
02791 }
02792 else{
02793 MPI_Send (buff,buffsize,MPI_LONG,0,myrank,MPI_COMM_WORLD);
02794 }
02795 MPI_Barrier (MPI_COMM_WORLD);
02796
02797 delete [] buff;
02798
02799 if (myrank == 0){
02800
02801
02802
02803
02804
02805
02806
02807
02808 a = gnn[0];
02809 gnn[0] = 0;
02810 for(i = 1; i < tnnp+1; i++){
02811 b = gnn[i];
02812 gnn[i] = gnn[i-1] + a;
02813 a = b;
02814 }
02815 tndof = gnn[tnnp];
02816
02817
02818 fprintf(out,"\n\n\ngnn array\n\n");
02819 for(i = 0; i < tnnp+1; i++){
02820 fprintf(out,"%ld %ld\n",i+1,gnn[i]);
02821 }
02822 }
02823 }
02824
02825
02826
02827
02828
02829
02830
02831
02832
02833
02834
02835
02836
02837
02838 void partop::assemble_pgcn(gtopology *top,long *domproc,FILE *out)
02839 {
02840 long i,ii,j,k,m;
02841 long adr,ndofn,buffsize;
02842 long *buff;
02843 MPI_Status stat;
02844
02845 buffsize = maxndof;
02846 buff = new long [buffsize];
02847
02848 m = 0;
02849 for (i = 0;i < nn; i++){
02850 ndofn = top->give_ndofn (i);
02851 for (j = 0; j < ndofn; j++){
02852 buff[m] = top -> give_dof (i,j);
02853 m++;
02854 }
02855 }
02856
02857
02858 if (myrank==0){
02859
02860 fprintf (stdout,"\n\n\n tndof je %ld\n\n\n",tndof);
02861
02862
02863
02864
02865
02866 if (pgcn == NULL){
02867 pgcn = new long[tndof];
02868 }
02869
02870
02871 for(i = 0; i < tndof; i++) pgcn[i] = -10;
02872
02873
02874 k = domproc[0];
02875 m = 0;
02876 for (i = 0;i < nnsd[k]; i++){
02877 adr = gnn[allnodes[k][i]];
02878
02879 ndofn = gnn[allnodes[k][i]+1] - adr;
02880
02881 for (j = 0;j < ndofn; j++){
02882 pgcn[adr+j] = buff[m];
02883 m++;
02884 }
02885 }
02886
02887
02888 for (i= 1 ;i < nproc; i++){
02889 MPI_Recv (buff,buffsize,MPI_LONG,MPI_ANY_SOURCE,MPI_ANY_TAG,MPI_COMM_WORLD,&stat);
02890
02891 k = domproc[stat.MPI_TAG];
02892 m = 0;
02893 for (ii = 0;ii< nnsd[k]; ii++){
02894 adr = gnn[allnodes[k][ii]];
02895 ndofn = gnn[allnodes[k][ii]+1] - adr;
02896 for (j = 0;j < ndofn; j++){
02897 pgcn[adr+j] = buff[m];
02898 m++;
02899 }
02900 }
02901 }
02902 }
02903 else{
02904 MPI_Send (buff,buffsize,MPI_LONG,0,myrank,MPI_COMM_WORLD);
02905 }
02906 MPI_Barrier (MPI_COMM_WORLD);
02907
02908 delete []buff;
02909
02910
02911
02912
02913
02914
02915
02916
02917
02918
02919
02920
02921
02922 if(myrank == 0){
02923 m = 1;
02924 for(i = 0; i < tndof; i++){
02925 if(pgcn[i] > 0){
02926 pgcn[i] = m;
02927 m++;
02928 }
02929 }
02930
02931 fprintf(out,"\n\n\n Global code numbers\n\n");
02932 for(i = 0; i < tndof; i++){
02933 fprintf(out,"%ld %ld\n",i,pgcn[i]);
02934 }
02935 tndof = m-1;
02936 fprintf(out,"\n\n\n number of DOFs on whole problem %ld\n\n",tndof);
02937 }
02938
02939
02940
02941 }
02942
02943
02944
02945
02946
02947
02948
02949
02950
02951
02952
02953
02954
02955
02956
02957
02958
02959
02960
02961
02962
02963 void partop::assemble_gcnd(gtopology *top,long *domproc,FILE *out)
02964 {
02965 long i,j,k,m;
02966 long ndofn,adr;
02967
02968 if(myrank == 0){
02969
02970 nud = new long[nproc];
02971 for(i = 0; i < nproc; i++){
02972 nud[i]=0;
02973 for(j = 0; j < nnsd[i]; j++){
02974 adr = gnn[allnodes[i][j]];
02975 ndofn = gnn[allnodes[i][j]+1] - adr;
02976 for (k = 0;k < ndofn; k++){
02977 if(pgcn[adr+k] > 0){
02978 nud[i]++;
02979 }
02980 }
02981 }
02982 }
02983
02984 fprintf(out,"\n\n\n Number of unknowns on subdomain\n");
02985 for(i = 0; i < nproc; i++){
02986 fprintf(out,"Domain %ld %ld\n",i+1,nud[i]);
02987 }
02988
02989 gcnd = new long*[nproc];
02990 for(i = 0;i < nproc; i++){
02991 gcnd[i] = new long[nud[i]];
02992 m = 0;
02993 for(j = 0; j < nnsd[i]; j++){
02994 adr = gnn[allnodes[i][j]];
02995 ndofn = gnn[allnodes[i][j]+1] - adr;
02996 for (k = 0;k < ndofn; k++){
02997 if(pgcn[adr+k] > 0){
02998 gcnd[i][m] = pgcn[adr+k] ;
02999 m++;
03000 }
03001 }
03002 }
03003 }
03004
03005
03006 fprintf(out,"\n\n\n Unknowns on subdomain gcnd\n");
03007 for(i = 0;i < nproc; i++){
03008 fprintf(out,"Domain %ld\n",i+1);
03009 for(j = 0; j < nud[i]; j++){
03010 fprintf(out,"%ld %ld\n",j,gcnd[i][j]);
03011 }
03012 }
03013
03014
03015 delete []gnn;
03016 delete []pgcn;
03017 }
03018
03019 }
03020
03021