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