00001 #include "gtopology.h"
00002 #include <stdio.h>
00003 #include <stdlib.h>
00004 #include <string.h>
00005 #include <limits.h>
00006 #include <math.h>
00007 #include "galias.h"
00008 #include "ordering.h"
00009 #include "gmatrix.h"
00010 #include "basefun.h"
00011 #include <time.h>
00012
00013 gtopology::gtopology ()
00014 {
00015
00016 nn=0;
00017
00018 ne=0;
00019
00020 nln=0;
00021
00022 ndof=0;
00023
00024 nidof=0;
00025
00026 nbdof=0;
00027
00028 nsad=0;
00029
00030
00031 nnso=0;
00032
00033
00034 neso=0;
00035
00036 ns=0;
00037
00038 nph=0;
00039
00040
00041 rst=0;
00042
00043 hangnod=0;
00044
00045
00046
00047 gnodes=NULL;
00048
00049 gelements=NULL;
00050
00051 lgnodes=NULL;
00052
00053 domsizes = NULL;
00054
00055
00056 nadjelnod=NULL;
00057
00058 adjelnod=NULL;
00059
00060 nadjelel=NULL;
00061
00062 adjelel=NULL;
00063
00064 nadjnodnod=NULL;
00065
00066 adjnodnod=NULL;
00067
00068
00069 bckcn=NULL; unln=NULL; unnl=NULL;
00070 cngtopcorr=NULL;
00071
00072 gphases=NULL;
00073
00074
00075
00076
00077
00078 dofcontr=0;
00079
00080
00081 cngen=1;
00082
00083 ordering=NULL;
00084
00085 nodren = no_renumbering;
00086
00087
00088
00089
00090
00091
00092 cnstate=0;
00093
00094 gf=NULL;
00095
00096 lnso=NULL;
00097
00098 leso=NULL;
00099
00100
00101
00102 edtype = noedge;
00103
00104 nged=0;
00105
00106 nser=0;
00107
00108 edgeser=NULL;
00109
00110 nedser=NULL;
00111
00112 edgelist = NULL;
00113
00114
00115 gedges=NULL;
00116
00117
00118 nen=0;
00119
00120 endnodes=NULL;
00121
00122
00123 stop=NULL;
00124
00125 nadjnodnod = NULL;
00126
00127 adjnodnod = NULL;
00128 }
00129
00130 gtopology::~gtopology ()
00131 {
00132 long i;
00133 delete [] gnodes;
00134 delete [] gelements;
00135 delete [] lgnodes;
00136
00137 delete [] domsizes;
00138
00139 delete [] nadjelnod;
00140 if (adjelnod)
00141 {
00142 for (i=0;i<nn;i++)
00143 delete [] adjelnod[i];
00144 }
00145 delete [] adjelnod;
00146
00147 delete [] nadjelel;
00148 if (adjelel)
00149 {
00150 for (i=0;i<ne;i++)
00151 delete [] adjelel[i];
00152 }
00153 delete [] adjelel;
00154
00155 delete [] bckcn;
00156 delete [] unln;
00157 delete [] unnl;
00158 delete [] cngtopcorr;
00159
00160 delete [] gphases;
00161 delete [] ordering;
00162
00163
00164
00165
00166 delete [] gf;
00167 delete [] lnso;
00168 delete [] leso;
00169
00170 delete [] nedser;
00171 delete [] edgeser;
00172 delete [] gedges;
00173
00174 if (edgelist)
00175 {
00176 for (i=0;i<nser;i++)
00177 delete [] edgelist[i];
00178 delete [] edgelist;
00179 }
00180
00181 delete [] endnodes;
00182
00183
00184 delete stop;
00185
00186 delete [] nadjnodnod;
00187
00188 if (adjnodnod)
00189 {
00190 for (i=0;i<nn;i++)
00191 delete [] adjnodnod[i];
00192 delete [] adjnodnod;
00193 }
00194 }
00195
00196
00197
00198
00199
00200
00201
00202
00203
00204
00205 void gtopology::initiate (long **icn,long *nc,long m)
00206 {
00207 long i;
00208 ne=m;
00209 cnstate=1;
00210 alloc_elements (ne);
00211 for (i=0;i<ne;i++){
00212 gelements[i].initiate (icn[i],nc[i]);
00213 }
00214 }
00215
00216
00217
00218
00219
00220
00221 void gtopology::alloc_nodes (long m)
00222 {
00223 nn=m;
00224 if (m<0){
00225 print_err("negative number of nodes is required", __FILE__, __LINE__, __func__);
00226 }
00227 gnodes = new gnode [nn];
00228 }
00229
00230
00231
00232
00233
00234
00235 void gtopology::alloc_elements (long m)
00236 {
00237 ne=m;
00238 if (m<0){
00239 print_err("negative number of elements is required", __FILE__, __LINE__, __func__);
00240 }
00241 gelements = new gelement [ne];
00242 }
00243
00244
00245
00246
00247
00248
00249 void gtopology::alloc_lnodes (long m)
00250 {
00251 nln=m;
00252 if (m<0){
00253 print_err("negative number of nodes is required", __FILE__, __LINE__, __func__);
00254 }
00255 lgnodes = new lgnode [nln];
00256 }
00257
00258
00259
00260
00261
00262
00263
00264
00265
00266 void gtopology::alloc_phases (long m)
00267 {
00268 nph=m;
00269 if (m<0){
00270 print_err("negative number of phases is required", __FILE__, __LINE__, __func__);
00271 }
00272 gphases = new gphase [nph];
00273 }
00274
00275
00276
00277
00278
00279
00280 long gtopology::give_ndofn (long nid)
00281 {
00282 long ndofn;
00283 ndofn = gnodes[nid].give_ndofn ();
00284 return ndofn;
00285 }
00286
00287
00288
00289
00290
00291
00292 long gtopology::give_nmult (long lnid)
00293 {
00294 return lgnodes[lnid].nmult;
00295 }
00296
00297
00298
00299
00300
00301
00302 long gtopology::give_ndofe (long eid)
00303 {
00304 return gelements[eid].give_ndofe ();
00305 }
00306
00307
00308
00309
00310
00311
00312
00313 long gtopology::give_gndofe (long eid)
00314 {
00315 long ndofe,nne;
00316 nne = gelements[eid].give_nne ();
00317 ndofe = gelements[eid].give_ndofe ();
00318 ndofe+= nne*2*gelements[eid].give_nmult ();
00319 return ndofe;
00320 }
00321
00322
00323
00324
00325
00326
00327 long gtopology::give_nnedge (long edid)
00328 {
00329 return gedges[edid].nn;
00330 }
00331
00332
00333
00334
00335
00336
00337 long gtopology::give_original_nne (long eid)
00338 {
00339 return gelements[eid].give_nne ();
00340 }
00341
00342
00343
00344
00345
00346
00347 long gtopology::give_nne (long eid)
00348 {
00349 long nne;
00350
00351 if (gelements[eid].give_nmne () == 0){
00352
00353 nne = gelements[eid].give_nne ();
00354 }else{
00355
00356 nne = gelements[eid].give_nmne ();
00357 }
00358 return nne;
00359 }
00360
00361
00362
00363
00364
00365
00366 long gtopology::give_nmne (long eid)
00367 {
00368 return gelements[eid].give_nmne ();
00369 }
00370
00371
00372
00373
00374
00375
00376 long gtopology::give_cne (long eid)
00377 {
00378 return gelements[eid].give_cne ();
00379 }
00380
00381
00382
00383
00384
00385
00386
00387 long gtopology::give_dof (long nid,long n)
00388 {
00389 return gnodes[nid].give_dof (n);
00390 }
00391
00392
00393
00394
00395
00396
00397
00398
00399 void gtopology::save_dof (long nid,long n,long num)
00400 {
00401 gnodes[nid].save_dof (n,num);
00402 }
00403
00404
00405
00406
00407
00408
00409
00410 void gtopology::give_nodes (long eid,ivector &nod)
00411 {
00412 if (gelements[eid].nmne==0){
00413
00414 gelements[eid].give_nodes (nod);
00415 }else{
00416
00417 gelements[eid].give_master_nodes (nod);
00418 }
00419 }
00420
00421
00422
00423
00424
00425
00426
00427 void gtopology::give_original_nodes (long eid,ivector &nod)
00428 {
00429 gelements[eid].give_nodes (nod);
00430 }
00431
00432
00433
00434
00435
00436
00437
00438
00439
00440 void gtopology::give_master_nodes (long eid,ivector &nod)
00441 {
00442 gelements[eid].give_master_nodes (nod);
00443 }
00444
00445
00446
00447
00448
00449
00450
00451
00452
00453
00454 void gtopology::give_node_coord1d (vector &x,long eid)
00455 {
00456 long i,nne;
00457 nne = gelements[eid].give_nne ();
00458 if (x.n!=nne){
00459 print_err("wrong size of array is required", __FILE__, __LINE__, __func__);
00460 }
00461 ivector nodes(nne);
00462 gelements[eid].give_nodes (nodes);
00463
00464 for (i=0;i<nne;i++){
00465 x[i]=gnodes[nodes[i]].x;
00466 }
00467 }
00468
00469
00470
00471
00472
00473
00474
00475
00476
00477
00478 void gtopology::give_node_coord2d (vector &x,vector &y,long eid)
00479 {
00480 long i,nne;
00481 nne = gelements[eid].give_nne ();
00482 if (x.n!=nne || y.n!=nne){
00483 print_err("wrong size of array is required", __FILE__, __LINE__, __func__);
00484 }
00485 ivector nodes(nne);
00486 gelements[eid].give_nodes (nodes);
00487
00488 for (i=0;i<nne;i++){
00489 x[i]=gnodes[nodes[i]].x;
00490 y[i]=gnodes[nodes[i]].y;
00491 }
00492 }
00493
00494
00495
00496
00497
00498
00499
00500
00501
00502 void gtopology::give_node_coord2dxz (vector &x,vector &z,long eid)
00503 {
00504 long i,nne;
00505 nne = gelements[eid].give_nne ();
00506 if (x.n!=nne || z.n!=nne){
00507 print_err("wrong size of array is required", __FILE__, __LINE__, __func__);
00508 }
00509 ivector nodes(nne);
00510 gelements[eid].give_nodes (nodes);
00511
00512 for (i=0;i<nne;i++){
00513 x[i]=gnodes[nodes[i]].x;
00514 z[i]=gnodes[nodes[i]].z;
00515 }
00516 }
00517
00518
00519
00520
00521
00522
00523
00524
00525
00526 void gtopology::give_node_coord3d (vector &x,vector &y,vector &z,long eid)
00527 {
00528 long i,nne;
00529 nne = gelements[eid].give_nne ();
00530 if (x.n!=nne || y.n!=nne || z.n!=nne){
00531 print_err("wrong size of array is assumed.",__FILE__,__LINE__,__func__);
00532 }
00533 ivector nodes(nne);
00534 gelements[eid].give_nodes (nodes);
00535
00536 for (i=0;i<nne;i++){
00537 x[i]=gnodes[nodes[i]].x;
00538 y[i]=gnodes[nodes[i]].y;
00539 z[i]=gnodes[nodes[i]].z;
00540 }
00541 }
00542
00543
00544
00545
00546
00547
00548
00549
00550
00551 void gtopology::give_code_numbers (long eid,long *cn)
00552 {
00553 long i,j,k,nne,ndofn;
00554
00555 if (cnstate==0){
00556 print_err("code numbers are required but not available",__FILE__,__LINE__,__func__);
00557 }
00558
00559 if (gelements[eid].cne==1){
00560
00561 for (i=0;i<give_ndofe (eid);i++){
00562 cn[i]=gelements[eid].cn[i];
00563 }
00564 }
00565 else{
00566
00567 k=0;
00568 nne=give_nne (eid);
00569 ivector nod(nne);
00570 give_nodes (eid,nod);
00571 for (i=0;i<nne;i++){
00572 ndofn=give_ndofn (nod[i]);
00573 for (j=0;j<ndofn;j++){
00574 cn[k]=give_dof (nod[i],j);
00575 k++;
00576 }
00577 }
00578 }
00579
00580 }
00581
00582
00583
00584
00585
00586
00587
00588
00589
00590 void gtopology::give_cn (long eid,long *cn)
00591 {
00592 long i;
00593
00594 if (cnstate==0){
00595 print_err("code numbers are required but not available",__FILE__,__LINE__,__func__);
00596 }
00597
00598 if (gelements[eid].cne==1 || gelements[eid].cne==2){
00599
00600 for (i=0;i<give_ndofe (eid);i++){
00601 cn[i]=gelements[eid].cn[i];
00602 }
00603 }
00604 }
00605
00606
00607
00608
00609
00610
00611
00612
00613
00614
00615 void gtopology::give_gcode_numbers (long eid,long *cn)
00616 {
00617 long i,j,k,nne,ndofe,ndofn,lnid,lid,nmult;
00618 ivector nod;
00619
00620 if (cnstate==0){
00621
00622 }
00623
00624 ndofe = give_ndofe (eid);
00625
00626 if (gelements[eid].cne==1){
00627
00628 for (i=0;i<ndofe;i++){
00629 cn[i]=gelements[eid].cn[i];
00630 }
00631 }
00632 else{
00633
00634 k=0;
00635 nne=give_nne (eid);
00636 allocv (nne,nod);
00637 give_nodes (eid,nod);
00638 for (i=0;i<nne;i++){
00639 ndofn=give_ndofn (nod[i]);
00640 for (j=0;j<ndofn;j++){
00641 cn[k]=give_dof (nod[i],j);
00642 k++;
00643 }
00644 }
00645 destrv (nod);
00646 }
00647
00648 if (gelements[eid].nmult>0){
00649 k=ndofe;
00650 nne=give_nne (eid);
00651 allocv (nne,nod);
00652 give_nodes (eid,nod);
00653 for (i=0;i<nne;i++){
00654 lnid=unln[nod[i]];
00655 lid=unnl[nod[i]];
00656 nmult=give_nmult (lnid);
00657 give_mult_code_numbers (lnid,lid,cn+k);
00658 k+=nmult;
00659 }
00660 for (i=0;i<nne;i++){
00661 lnid=unln[nod[i]];
00662 lid=unnl[nod[i]]+1;
00663 nmult=give_nmult (lnid);
00664 give_mult_code_numbers (lnid,lid,cn+k);
00665 k+=nmult;
00666 }
00667 destrv (nod);
00668 }
00669
00670 }
00671
00672
00673
00674
00675
00676
00677
00678
00679
00680 void gtopology::give_node_code_numbers (long nid,long *cn)
00681 {
00682 long i,ndofn;
00683
00684 if (cnstate==0){
00685 print_err("code numbers are required but not available", __FILE__, __LINE__, __func__);
00686 }
00687
00688 ndofn = give_ndofn (nid);
00689 for (i=0;i<ndofn;i++){
00690 cn[i]=gnodes[nid].cn[i];
00691 }
00692 }
00693
00694
00695
00696
00697
00698
00699
00700
00701
00702
00703 void gtopology::give_mult_code_numbers (long nid,long lid,long *cn)
00704 {
00705 long i,nmult;
00706
00707 if (cnstate==0){
00708 print_err("code numbers are required but not available", __FILE__, __LINE__, __func__);
00709 }
00710
00711 nmult = give_nmult (nid);
00712 for (i=0;i<nmult;i++){
00713 cn[i]=lgnodes[nid].cn[lid][i];
00714 }
00715 }
00716
00717
00718
00719
00720
00721
00722
00723
00724
00725
00726
00727 void gtopology::give_endnode_code_numbers (long enid,long *ncn1,long *ncn2,long *mcn)
00728 {
00729 long ndofn,*nid;
00730
00731 nid = new long [2];
00732
00733
00734 endnodes[enid].give_node_numbers (nid);
00735
00736
00737 ndofn = give_ndofn (nid[0]);
00738
00739
00740 give_node_code_numbers (nid[0],ncn1);
00741 give_node_code_numbers (nid[1],ncn2);
00742
00743
00744 endnodes[enid].give_mult_code_numbers (mcn);
00745
00746 delete [] nid;
00747 }
00748
00749
00750
00751
00752
00753
00754
00755
00756
00757
00758
00759
00760
00761
00762 void gtopology::give_edge_code_numbers (long edid,long fln,long *ncn1,long *ncn2,long *mcn)
00763 {
00764 long ndofn,*nid;
00765
00766 nid = new long [2];
00767
00768 if (fln==1){
00769 gedges[edid].give_first_node_numbers (nid);
00770 }
00771 if (fln==2){
00772 gedges[edid].give_last_node_numbers (nid);
00773 }
00774
00775
00776 ndofn = give_ndofn (nid[0]);
00777
00778
00779 give_node_code_numbers (nid[0],ncn1);
00780 give_node_code_numbers (nid[1],ncn2);
00781
00782 gedges[edid].give_mult_code_numbers (fln,mcn);
00783
00784 delete [] nid;
00785 }
00786
00787
00788
00789
00790
00791
00792
00793
00794
00795
00796 gelemtype gtopology::give_elem_type (long eid)
00797 {
00798 return gelements[eid].get;
00799 }
00800
00801
00802
00803
00804
00805
00806
00807
00808
00809 void gtopology::give_end_nodes (long eid,long *nodes)
00810 {
00811 long nne;
00812 gelemtype get;
00813 ivector nod;
00814 long *enod;
00815
00816 enod=new long [2];
00817
00818
00819 nne = give_original_nne (eid);
00820
00821 allocv (nne,nod);
00822
00823 give_original_nodes (eid,nod);
00824
00825 get = give_elem_type (eid);
00826
00827 switch (get){
00828 case linbar:{
00829 linbar_endpoints (enod);
00830 nodes[0]=nod[enod[0]];
00831 nodes[1]=nod[enod[1]];
00832 break;
00833 }
00834 case quadbar:{
00835 quadbar_endpoints (enod);
00836 nodes[0]=nod[enod[0]];
00837 nodes[1]=nod[enod[1]];
00838 break;
00839 }
00840
00841 default:{
00842 print_err("unknown element type is required", __FILE__, __LINE__, __func__);
00843 }
00844 }
00845
00846 delete [] enod;
00847 destrv (nod);
00848 }
00849
00850
00851
00852
00853
00854
00855
00856
00857
00858
00859
00860
00861
00862
00863 void gtopology::give_edge_loc_nodes (long eid,long edid,ivector &edgenod)
00864 {
00865 gelemtype get;
00866
00867
00868 get = give_elem_type (eid);
00869
00870 switch (get)
00871 {
00872 case linbar:
00873 edgenod[0] = 0;
00874 edgenod[1] = 1;
00875 break;
00876 case quadbar:
00877 edgenod[0] = 0;
00878 edgenod[1] = 1;
00879 edgenod[2] = 2;
00880 break;
00881 case lintriag:
00882 lintriangle_edgnod(edgenod.a,edid);
00883 break;
00884 case quadtriag:
00885 quadtriangle_edgnod(edgenod.a,edid);
00886 break;
00887 case linquad:
00888 linquadrilat_edgnod (edgenod.a,edid);
00889 break;
00890 case quadquad:
00891 quadquadrilat_edgnod (edgenod.a,edid);
00892 break;
00893 case lintetra:
00894 lintetrahedral_edgnod(edgenod.a,edid);
00895 break;
00896 case quadtetra:
00897 quadtetrahedral_edgnod (edgenod.a,edid);
00898 break;
00899 case linhexa:
00900 linhexahedral_edgnod(edgenod.a,edid);
00901 break;
00902 case quadhexa:
00903 quadhexahedral_edgnod (edgenod.a,edid);
00904 break;
00905 default:
00906 print_err("unknown element type is required", __FILE__, __LINE__, __func__);
00907 }
00908 }
00909
00910
00911
00912
00913
00914
00915
00916
00917
00918
00919
00920
00921 void gtopology::give_edge_nodes (long eid,long edid,long *nodes)
00922 {
00923 long i,nne,nned;
00924 gelemtype get;
00925 ivector nod;
00926 long *edgenod;
00927
00928
00929 nne = give_original_nne (eid);
00930
00931 allocv (nne,nod);
00932
00933 give_original_nodes (eid,nod);
00934
00935 get = give_elem_type (eid);
00936
00937 nned = give_nned (eid);
00938
00939 edgenod = new long[nned];
00940 switch (get){
00941 case lintriag:{
00942 lintriangle_edgnod(edgenod,edid);
00943 for(i = 0; i < nned; i++){
00944 nodes[i]=nod[edgenod[i]];
00945 }
00946 break;
00947 }
00948 case quadtriag:{
00949 quadtriangle_edgnod(edgenod,edid);
00950 for(i = 0; i < nned; i++){
00951 nodes[i]=nod[edgenod[i]];
00952 }
00953 break;
00954 }
00955 case linquad:{
00956 linquadrilat_edgnod (edgenod,edid);
00957 for(i = 0; i < nned; i++){
00958 nodes[i]=nod[edgenod[i]];
00959 }
00960 break;
00961 }
00962 case quadquad:{
00963 quadquadrilat_edgnod (edgenod,edid);
00964 for(i = 0; i < nned; i++){
00965 nodes[i]=nod[edgenod[i]];
00966 }
00967 break;
00968 }
00969 case lintetra:{
00970 lintetrahedral_edgnod(edgenod,edid);
00971 for(i = 0; i < nned; i++){
00972 nodes[i]=nod[edgenod[i]];
00973 }
00974 break;
00975 }
00976 case quadtetra:{
00977 quadtetrahedral_edgnod (edgenod,edid);
00978 for(i = 0; i < nned; i++){
00979 nodes[i]=nod[edgenod[i]];
00980 }
00981 break;
00982 }
00983 case linhexa:{
00984 linhexahedral_edgnod(edgenod,edid);
00985 for(i = 0; i < nned; i++){
00986 nodes[i]=nod[edgenod[i]];
00987 }
00988 break;
00989 }
00990 case quadhexa:{
00991 quadhexahedral_edgnod (edgenod,edid);
00992 for(i = 0; i < nned; i++){
00993 nodes[i]=nod[edgenod[i]];
00994 }
00995 break;
00996 }
00997 default:{
00998 print_err("unknown element type is required", __FILE__, __LINE__, __func__);
00999 }
01000 }
01001
01002 destrv (nod);
01003 delete []edgenod;
01004 }
01005
01006
01007
01008
01009
01010
01011
01012
01013
01014
01015 void gtopology::give_surf_nodes (long eid,long surfid,long *nodes)
01016 {
01017 long i,nne,nnsurf;
01018 gelemtype get;
01019 ivector nod;
01020 long *surfnod;
01021
01022
01023 nne = give_nne (eid);
01024
01025 allocv (nne,nod);
01026
01027 give_nodes (eid,nod);
01028
01029 get = give_elem_type (eid);
01030
01031 nnsurf = give_nnsurf (eid);
01032
01033 surfnod = new long[nnsurf];
01034 switch(get){
01035 case lintetra:{
01036 lintetrahedral_surfnod (surfnod,surfid);
01037 for(i = 0; i < nnsurf; i++){
01038 nodes[i]=nod[surfnod[i]];
01039 }
01040 break;
01041 }
01042 case quadtetra:{
01043 quadtetrahedral_surfnod (surfnod,surfid);
01044 for(i = 0; i < nnsurf; i++){
01045 nodes[i]=nod[surfnod[i]];
01046 }
01047 break;
01048 }
01049 case linhexa:{
01050 linhexahedral_surfnod (surfnod,surfid);
01051 for(i = 0; i < nnsurf; i++){
01052 nodes[i]=nod[surfnod[i]];
01053 }
01054 break;
01055 }
01056 case quadhexa:{
01057 quadhexahedral_surfnod (surfnod,surfid);
01058 for(i = 0; i < nnsurf; i++){
01059 nodes[i]=nod[surfnod[i]];
01060 }
01061 break;
01062 }
01063
01064 default:{
01065 print_err("unknown element type is required", __FILE__, __LINE__, __func__);
01066 }
01067 }
01068 destrv (nod);
01069 delete []surfnod;
01070 }
01071
01072
01073
01074
01075
01076
01077
01078
01079 long gtopology::give_nen (long eid)
01080 {
01081 long ned=0;
01082 gelemtype get;
01083
01084 get = give_elem_type (eid);
01085
01086 switch (get){
01087 case linbar:{ ned=2; break; }
01088 case quadbar:{ ned=2; break; }
01089 case lintriag:{ ned=0; break; }
01090 case quadtriag:{ ned=0; break; }
01091 case linquad:{ ned=0; break; }
01092 case quadquad:{ ned=0; break; }
01093 case lintetra:{ ned=0; break; }
01094 case quadtetra:{ ned=0; break; }
01095 case linhexa:{ ned=0; break; }
01096 case quadhexa:{ ned=0; break; }
01097 default:{
01098 print_err("unknown element type is required", __FILE__, __LINE__, __func__);
01099 }
01100 }
01101 return ned;
01102 }
01103
01104
01105
01106
01107
01108
01109
01110
01111 long gtopology::give_ned (long eid)
01112 {
01113 long ned=0;
01114 gelemtype get;
01115
01116 get = give_elem_type (eid);
01117
01118 switch (get){
01119 case linbar:{ ned=0; break; }
01120 case quadbar:{ ned=0; break; }
01121 case lintriag:{ ned=3; break; }
01122 case quadtriag:{ ned=3; break; }
01123 case linquad:{ ned=4; break; }
01124 case quadquad:{ ned=4; break; }
01125 case lintetra:{ ned=6; break; }
01126 case quadtetra:{ ned=6; break; }
01127 case linhexa:{ ned=12; break; }
01128 case quadhexa:{ ned=12; break; }
01129 default:{
01130 print_err("unknown element type is required", __FILE__, __LINE__, __func__);
01131 }
01132 }
01133 return ned;
01134 }
01135
01136
01137
01138
01139
01140
01141
01142
01143 long gtopology::give_nned (long eid)
01144 {
01145 long nned=0;
01146 gelemtype get;
01147
01148 get = give_elem_type (eid);
01149
01150 switch (get){
01151 case linbar:{ nned=2; break; }
01152 case quadbar:{ nned=3; break; }
01153 case lintriag:{ nned=2; break; }
01154 case quadtriag:{ nned=3; break; }
01155 case linquad:{ nned=2; break; }
01156 case quadquad:{ nned=3; break; }
01157 case lintetra:{ nned=2; break; }
01158 case quadtetra:{ nned=3; break; }
01159 case linhexa:{ nned=2; break; }
01160 case quadhexa:{ nned=3; break; }
01161 default:{
01162 print_err("unknown element type is required", __FILE__, __LINE__, __func__);
01163 }
01164 }
01165 return nned;
01166 }
01167
01168
01169
01170
01171
01172
01173
01174
01175 long gtopology::give_nsurf (long eid)
01176 {
01177 long nsurf;
01178 gelemtype get;
01179
01180 get = give_elem_type (eid);
01181
01182 switch (get){
01183 case lintetra:{ nsurf=4; break; }
01184 case quadtetra:{ nsurf=4; break; }
01185 case linhexa:{ nsurf=6; break; }
01186 case quadhexa:{ nsurf=6; break; }
01187 default:{
01188 print_err("unknown element type is required", __FILE__, __LINE__, __func__);
01189 }
01190 }
01191 return nsurf;
01192 }
01193
01194
01195
01196
01197
01198
01199
01200
01201 long gtopology::give_nnsurf (long eid)
01202 {
01203 long nnsurf;
01204 gelemtype get;
01205
01206 get = give_elem_type (eid);
01207
01208 switch (get){
01209 case lintetra:{ nnsurf=3; break; }
01210 case linhexa:{ nnsurf=4; break; }
01211 case quadtetra:{ nnsurf=6; break; }
01212 case quadhexa:{ nnsurf=8; break; }
01213 default:{
01214 print_err("unknown element type is required", __FILE__, __LINE__, __func__);
01215 }
01216 }
01217 return nnsurf;
01218 }
01219
01220
01221
01222
01223
01224
01225
01226
01227
01228 long gtopology::give_degree (long eid)
01229 {
01230 long deg=0;
01231 gelemtype get;
01232
01233 get = give_elem_type (eid);
01234
01235 switch (get){
01236 case linbar:{ deg=1; break; }
01237 case quadbar:{ deg=2; break; }
01238 case lintriag:{ deg=1; break; }
01239 case quadtriag:{ deg=2; break; }
01240 case linquad:{ deg=1; break; }
01241 case quadquad:{ deg=2; break; }
01242 case lintetra:{ deg=1; break; }
01243 case quadtetra:{ deg=2; break; }
01244 case linhexa:{ deg=1; break; }
01245 case quadhexa:{ deg=2; break; }
01246 default:{
01247 print_err("unknown element type is required", __FILE__, __LINE__, __func__);
01248 }
01249 }
01250 return deg;
01251 }
01252
01253
01254
01255
01256
01257
01258 long gtopology::give_whole_dim(long eid)
01259 {
01260
01261 gelemtype get;
01262 long dim;
01263 get = give_elem_type (eid);
01264 switch(get){
01265 case linbar:{ dim = 1; break; }
01266 case quadbar:{ dim = 1; break; }
01267 case lintriag:{ dim = 2; break; }
01268 case quadtriag:{ dim = 2; break; }
01269 case linquad:{ dim = 2; break; }
01270 case quadquad:{ dim = 2; break; }
01271 case lintetra:{ dim = 3; break; }
01272 case quadtetra:{ dim = 3; break; }
01273 case linhexa:{ dim = 3; break; }
01274 case quadhexa:{ dim = 3; break; }
01275 default:{
01276 print_err("unknown element type is required\n", __FILE__, __LINE__, __func__);
01277 }
01278 }
01279 return (dim);
01280 }
01281
01282
01283
01284
01285
01286
01287
01288
01289
01290
01291 long gtopology::codenum_generation (FILE *out)
01292 {
01293 long i,ndof;
01294
01295 ordering = new long [nn];
01296
01297
01298 switch (nodren){
01299 case no_renumbering:{
01300 for (i=0;i<nn;i++){
01301 ordering[i]=i;
01302 }
01303 break;
01304 }
01305 case cuthill_mckee:
01306 case rev_cuthill_mckee:{
01307 cuthill_mckee_renumb (out);
01308 break;
01309 }
01310 case sloan:{
01311 sloan_renumb(out);
01312 break;
01313 }
01314 default:{
01315 print_err("unknown type of renumbering is required",__FILE__,__LINE__,__func__);
01316 }
01317 }
01318
01319
01320 if (cngen==0){
01321 ndof = 0;
01322 }
01323
01324 if (cngen==1){
01325
01326 ndof = gencodnum ();
01327 }
01328 if (cngen==2){
01329
01330
01331
01332
01333 ndof = schur_ordering ();
01334 }
01335 if (cngen==3){
01336
01337
01338
01339 ndof = saddlepoint_ordering (stop->ns,stop->nnsd,stop->ltg);
01340 }
01341
01342 delete [] ordering;
01343 ordering = NULL;
01344
01345 if (hangnod==1){
01346
01347
01348 dof_transfer_hangnod (out);
01349 }
01350
01351 return ndof;
01352 }
01353
01354
01355
01356
01357
01358
01359 long gtopology::gencodnum ()
01360 {
01361 long i,j,k,ii,ndofn,*aux;
01362
01363
01364 ndof=0;
01365 for (i=0;i<nn;i++){
01366 ndofn=give_ndofn (i);
01367 for (j=0;j<ndofn;j++){
01368 k=give_dof (i,j);
01369 if (k>ndof)
01370 ndof=k;
01371 }
01372 }
01373 ndof--;
01374 if (ndof<0) ndof=0;
01375 aux = new long [ndof];
01376 for (i=0;i<ndof;i++){
01377 aux[i]=-1;
01378 }
01379
01380
01381 ndof=1;
01382 for (i=0;i<nn;i++){
01383 ii=ordering[i];
01384 ndofn=give_ndofn (ii);
01385 for (j=0;j<ndofn;j++){
01386 k=give_dof (ii,j);
01387 if (k<0) continue;
01388 if (k==0) continue;
01389 if (k==1){
01390 gnodes[ii].cn[j]=ndof; ndof++;
01391 }
01392 if (k>1){
01393 if (aux[k-2]==-1){
01394 gnodes[ii].cn[j]=ndof;
01395 aux[k-2]=ndof;
01396 ndof++;
01397 }
01398 else{
01399 gnodes[ii].cn[j]=aux[k-2];
01400 }
01401 }
01402 }
01403 }
01404 ndof--;
01405
01406 delete [] aux;
01407
01408
01409 for (i=0;i<ne;i++){
01410 if (gelements[i].cne==1){
01411 for (j=0;j<give_ndofe (i);j++){
01412 if (gelements[i].cn[j]>ndof)
01413 ndof=gelements[i].cn[j];
01414 }
01415 }
01416 }
01417
01418
01419 cnstate=1;
01420
01421 return ndof;
01422 }
01423
01424
01425
01426
01427
01428
01429
01430
01431
01432
01433
01434 long gtopology::schur_ordering ()
01435 {
01436 long i,j,k,l,ii,jj,kk,ndofn,*aux;
01437
01438
01439 ndof=0;
01440 for (i=0;i<nn;i++){
01441 ndofn=give_ndofn (i);
01442 for (j=0;j<ndofn;j++){
01443 k=give_dof (i,j);
01444 if (k>ndof) ndof=k;
01445 }
01446 }
01447 ndof--;
01448 if (ndof<0) ndof=0;
01449 aux = new long [ndof];
01450 for (i=0;i<ndof;i++){
01451 aux[i]=-1;
01452 }
01453
01454
01455 jj=-1; kk=-1;
01456 ndof=1;
01457 for (l=0;l<ns;l++){
01458 for (i=0;i<stop->nnsd[l];i++){
01459 jj++;
01460 if (stop->ltg[l][i]==-1){
01461 ii=ordering[jj];
01462 ndofn=give_ndofn (ii);
01463 for (j=0;j<ndofn;j++){
01464 k=give_dof (ii,j);
01465 if (k<0) continue;
01466 if (k==0) continue;
01467 if (k==1){
01468 gnodes[ii].cn[j]=ndof; ndof++;
01469 }
01470 if (k>1){
01471 if (aux[k-2]==-1){
01472 gnodes[ii].cn[j]=ndof;
01473 aux[k-2]=ndof;
01474 ndof++;
01475 }
01476 else{
01477 gnodes[ii].cn[j]=aux[k-2];
01478 }
01479 }
01480 }
01481 }
01482 }
01483
01484
01485 nidof=ndof-1;
01486
01487
01488 jj=kk;
01489 for (i=0;i<stop->nnsd[l];i++){
01490 jj++;
01491 if (stop->ltg[l][i]>-1){
01492 ii=ordering[jj];
01493 ndofn=give_ndofn (ii);
01494 for (j=0;j<ndofn;j++){
01495 k=give_dof (ii,j);
01496 if (k<0) continue;
01497 if (k==0) continue;
01498 if (k==1){
01499 gnodes[ii].cn[j]=ndof; ndof++;
01500 }
01501 if (k>1){
01502 if (aux[k-2]==-1){
01503 gnodes[ii].cn[j]=ndof;
01504 aux[k-2]=ndof;
01505 ndof++;
01506 }
01507 else{
01508 gnodes[ii].cn[j]=aux[k-2];
01509 }
01510 }
01511 }
01512 }
01513 }
01514 kk+=stop->nnsd[l];
01515 }
01516 ndof--;
01517
01518 delete [] aux;
01519
01520
01521 for (i=0;i<ne;i++){
01522 if (gelements[i].cne==1){
01523 for (j=0;j<give_ndofe (i);j++){
01524 if (gelements[i].cn[j]>ndof)
01525 ndof=gelements[i].cn[j];
01526 }
01527 }
01528 }
01529
01530
01531 nbdof=ndof-nidof;
01532
01533
01534 cnstate=1;
01535
01536 return ndof;
01537
01538
01539
01540
01541
01542
01543
01544
01545
01546
01547
01548
01549
01550
01551
01552
01553
01554
01555
01556
01557
01558
01559
01560
01561
01562
01563
01564
01565
01566
01567
01568
01569
01570
01571
01572
01573
01574
01575
01576
01577
01578
01579
01580
01581
01582
01583
01584
01585
01586
01587
01588
01589
01590
01591
01592
01593
01594
01595
01596
01597
01598
01599
01600
01601
01602
01603
01604
01605
01606
01607
01608
01609
01610
01611
01612
01613
01614
01615
01616
01617
01618
01619
01620
01621
01622
01623
01624
01625
01626
01627
01628
01629
01630
01631
01632
01633
01634 }
01635
01636
01637
01638
01639
01640
01641
01642
01643
01644
01645
01646
01647
01648
01649
01650
01651
01652
01653 long gtopology::saddlepoint_ordering (long ns,long *nnsd,long **ltg)
01654 {
01655 long i,j,k,ii,jj,nid,ndofn,*aux;
01656
01657
01658 ndof=0;
01659 for (i=0;i<nn;i++){
01660 ndofn=give_ndofn (i);
01661 for (j=0;j<ndofn;j++){
01662 k=give_dof (i,j);
01663 if (k>ndof) ndof=k;
01664 }
01665 }
01666 ndof--;
01667 if (ndof<0) ndof=0;
01668 aux = new long [ndof];
01669 for (i=0;i<ndof;i++){
01670 aux[i]=-1;
01671 }
01672
01673
01674 nid=0;
01675 ndof=1;
01676 for (jj=0;jj<ns;jj++){
01677 for (i=0;i<nnsd[jj];i++){
01678 if (ltg[jj][i]==-1){
01679 ii=ordering[nid];
01680 ndofn=give_ndofn (ii);
01681 for (j=0;j<ndofn;j++){
01682 k=give_dof (ii,j);
01683 if (k<0) continue;
01684 if (k==0) continue;
01685 if (k==1){
01686 gnodes[ii].cn[j]=ndof; ndof++;
01687 }
01688 if (k>1){
01689 if (aux[k-2]==-1){
01690 gnodes[ii].cn[j]=ndof;
01691 aux[k-2]=ndof;
01692 ndof++;
01693 }
01694 else{
01695 gnodes[ii].cn[j]=aux[k-2];
01696 }
01697 }
01698 }
01699 }
01700 nid++;
01701 }
01702 }
01703
01704 nidof=ndof-1;
01705
01706
01707 nid=0;
01708 for (jj=0;jj<ns;jj++){
01709 for (i=0;i<nnsd[jj];i++){
01710 if (ltg[jj][i]!=-1){
01711 ii=ordering[nid];
01712 ndofn=give_ndofn (ii);
01713 for (j=0;j<ndofn;j++){
01714 k=give_dof (ii,j);
01715 if (k<0) continue;
01716 if (k==0) continue;
01717 if (k==1){
01718 gnodes[ii].cn[j]=ndof; ndof++;
01719 }
01720 if (k>1){
01721 if (aux[k-2]==-1){
01722 gnodes[ii].cn[j]=ndof;
01723 aux[k-2]=ndof;
01724 ndof++;
01725 }
01726 else{
01727 gnodes[ii].cn[j]=aux[k-2];
01728 }
01729 }
01730 }
01731 }
01732 nid++;
01733 }
01734 }
01735
01736
01737 delete [] aux;
01738
01739
01740 for (i=0;i<ne;i++){
01741 if (gelements[i].cne==1){
01742 for (j=0;j<give_ndofe (i);j++){
01743 if (gelements[i].cn[j]>ndof)
01744 ndof=gelements[i].cn[j];
01745 }
01746 }
01747 }
01748
01749
01750 nbdof=ndof-nidof;
01751
01752
01753
01754
01755 ndof = codenum_multip ();
01756
01757 ndof--;
01758
01759
01760 cnstate=1;
01761
01762 fprintf (stdout,"\n number of DOFs %ld",ndof);
01763 fprintf (stdout,"\n number of internal DOFs %ld",nidof);
01764 fprintf (stdout,"\n number of interface DOFs %ld",nbdof);
01765 fprintf (stdout,"\n number of unknowns in sadd. point. prob. %ld",nsad);
01766
01767
01768 return ndof;
01769 }
01770
01771
01772
01773
01774
01775
01776
01777
01778
01779
01780
01781
01782
01783
01784
01785
01786
01787
01788
01789
01790
01791
01792
01793
01794
01795 void gtopology::comptop (gtopology *gt1,gtopology *gt2)
01796 {
01797 long i,j,k,l,nn1,nn2,ne1,ne2,ndofn1,nne;
01798 ivector nod;
01799
01800
01801 nn1 = gt1->nn;
01802 nn2 = gt2->nn;
01803 if (nn1!=nn2){
01804 print_err("two topologies have different numbers of nodes", __FILE__, __LINE__, __func__);
01805 abort ();
01806 }
01807
01808 nn = nn1;
01809 gnodes = new gnode [nn];
01810
01811
01812
01813 for (i=0;i<nn;i++){
01814 ndofn1 = gt1->gnodes[i].ndofn;
01815 gnodes[i].ndofn = ndofn1;
01816
01817 gnodes[i].x=gt1->gnodes[i].x;
01818 gnodes[i].y=gt1->gnodes[i].y;
01819 gnodes[i].z=gt1->gnodes[i].z;
01820 }
01821
01822
01823 ne1 = gt1->ne;
01824 ne2 = gt2->ne;
01825 if (ne1!=ne2){
01826 print_err("two topologies have different numbers of elements", __FILE__, __LINE__, __func__);
01827 abort ();
01828 }
01829
01830 ne = ne1;
01831 gelements = new gelement [ne];
01832
01833
01834 for (i=0;i<ne;i++){
01835 nne = gt2->give_nne (i);
01836 allocv (nne,nod);
01837 gt2->give_nodes (i,nod);
01838 for (j=0;j<nne;j++){
01839 if (gt2->gnodes[nod[j]].ai==1) continue;
01840 gnodes[nod[j]].ndofn+=gt2->gnodes[nod[j]].ndofn;
01841 gt2->gnodes[nod[j]].ai=1;
01842 }
01843 destrv (nod);
01844 }
01845
01846 for (i=0;i<nn;i++){
01847 gt2->gnodes[i].ai=0;
01848 }
01849
01850
01851 for (i=0;i<nn;i++){
01852 gnodes[i].cn = new long [gnodes[i].ndofn];
01853 for (j=0;j<gt1->gnodes[i].ndofn;j++){
01854 gnodes[i].cn[j]=gt1->gnodes[i].cn[j];
01855 }
01856 }
01857
01858
01859 for (i=0;i<ne;i++){
01860 nne = gt2->give_nne (i);
01861 allocv (nne,nod);
01862 gt2->give_nodes (i,nod);
01863 for (j=0;j<nne;j++){
01864 if (gt2->gnodes[nod[j]].ai==1) continue;
01865 l=gt1->gnodes[nod[j]].ndofn;
01866 for (k=0;k<gt2->gnodes[nod[j]].ndofn;k++){
01867 gnodes[nod[j]].cn[l]=gt2->gnodes[nod[j]].cn[k]; l++;
01868 }
01869 gt2->gnodes[nod[j]].ai=1;
01870 }
01871 destrv (nod);
01872 }
01873
01874
01875
01876 for (i=0;i<ne1;i++){
01877 gelements[i].ndofe = gt1->gelements[i].ndofe + gt2->gelements[i].ndofe;
01878 nne = gt1->gelements[i].nne;
01879 gelements[i].nne = nne;
01880 gelements[i].nodes = new long [nne];
01881 for (j=0;j<nne;j++){
01882 gelements[i].nodes[j] = gt1->gelements[i].nodes[j];
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 void gtopology::cndistr (gtopology *gt1,gtopology *gt2)
01911 {
01912 long i,j,k,ndofn1,ndofn2;
01913
01914 for (i=0;i<nn;i++){
01915 ndofn1=gt1->gnodes[i].ndofn;
01916 if (gnodes[i].ndofn>ndofn1){
01917 for (j=0;j<ndofn1;j++){
01918 gt1->gnodes[i].cn[j]=gnodes[i].cn[j];
01919 }
01920 k=ndofn1;
01921 ndofn2=gt2->gnodes[i].ndofn;
01922 for (j=0;j<ndofn2;j++){
01923 gt2->gnodes[i].cn[j]=gnodes[i].cn[k];
01924 k++;
01925 }
01926 }
01927 else{
01928 for (j=0;j<ndofn1;j++){
01929 gt1->gnodes[i].cn[j]=gnodes[i].cn[j];
01930 }
01931 }
01932 }
01933
01934 }
01935
01936
01937
01938
01939
01940
01941
01942
01943
01944 void gtopology::adjacelem (FILE *out)
01945 {
01946 long i,j,k,min,ind,s,nne,last;
01947 long **aux;
01948 ivector nodes;
01949
01950
01951 if (nadjelnod==NULL)
01952 nadjelnod = new long [nn];
01953 memset (nadjelnod,0,nn*sizeof(long));
01954
01955
01956 for (i=0;i<ne;i++){
01957 nne = give_nne (i);
01958 reallocv (nne,nodes);
01959 give_nodes (i,nodes);
01960 for (j=0;j<nne;j++){
01961 nadjelnod[nodes[j]]++;
01962 }
01963 }
01964
01965
01966 if (adjelnod==NULL){
01967 adjelnod = new long* [nn];
01968 for (i=0;i<nn;i++){
01969 adjelnod[i] = new long [nadjelnod[i]];
01970 }
01971 }
01972 for (i=0;i<nn;i++){
01973 nadjelnod[i] = 0;
01974 }
01975
01976
01977
01978 for (i=0;i<ne;i++){
01979 nne = give_nne (i);
01980 reallocv (nne,nodes);
01981 give_nodes (i,nodes);
01982 for (j=0;j<nne;j++){
01983 adjelnod[nodes[j]][nadjelnod[nodes[j]]++] = i;
01984 }
01985 }
01986
01987 if (out)
01988 {
01989
01990 fprintf (out,"\n\n\n NUMBERS OF ELEMENTS ADJACENT TO NODES\n");
01991 for (i=0;i<nn;i++){
01992 fprintf (out,"\n node %5ld %ld",i,nadjelnod[i]);
01993 }
01994
01995 fprintf (out,"\n\n\n LIST OF ELEMENTS ADJACENT TO NODES\n");
01996 for (i=0;i<nn;i++){
01997 fprintf (out,"\n node %4ld ",i);
01998 for (j=0;j<nadjelnod[i];j++){
01999 fprintf (out," %ld",adjelnod[i][j]);
02000 }
02001 }
02002 fprintf (out,"\n\n\n");
02003 }
02004
02005
02006
02007 nadjelel = new long [ne];
02008 memset (nadjelel,0,ne*sizeof(long));
02009
02010
02011 for (i=0;i<ne;i++){
02012 nne=give_nne (i);
02013 reallocv (nne,nodes);
02014 give_nodes (i,nodes);
02015 for (j=0;j<nne;j++){
02016 nadjelel[i]+=nadjelnod[nodes[j]];
02017 }
02018 }
02019
02020
02021 aux = new long* [ne];
02022 for (i=0;i<ne;i++){
02023 aux[i] = new long [nadjelel[i]];
02024 }
02025
02026
02027
02028 memset (nadjelel,0,ne*sizeof(long));
02029 for (i=0;i<ne;i++){
02030 nne=give_nne (i);
02031 reallocv (nne,nodes);
02032 give_nodes (i,nodes);
02033 for (j=0;j<nne;j++){
02034 for (k=0;k<nadjelnod[nodes[j]];k++){
02035 aux[i][nadjelel[i]]=adjelnod[nodes[j]][k];
02036 nadjelel[i]++;
02037 }
02038 }
02039 }
02040
02041
02042 for (i=0;i<ne;i++){
02043 last=-1;
02044 for (j=0;j<nadjelel[i];j++){
02045 min=LONG_MAX;
02046 for (k=j;k<nadjelel[i];k++){
02047 if (min>aux[i][k]){
02048 min=aux[i][k]; ind=k;
02049 }
02050 }
02051 if (last==min){
02052 aux[i][ind]=aux[i][nadjelel[i]-1];
02053 nadjelel[i]--; j--;
02054 }
02055 else{
02056 s=aux[i][j];
02057 aux[i][j]=min;
02058 aux[i][ind]=s;
02059 last=min;
02060 }
02061 }
02062 }
02063
02064
02065 adjelel = new long* [ne];
02066 for (i=0;i<ne;i++){
02067 adjelel[i] = new long [nadjelel[i]];
02068 memset (adjelel[i],0,nadjelel[i]*sizeof(long));
02069 }
02070
02071 for (i=0;i<ne;i++){
02072 for (j=0;j<nadjelel[i];j++){
02073 adjelel[i][j]=aux[i][j];
02074 }
02075 }
02076
02077 if (out)
02078 {
02079
02080 fprintf (out,"\n\n\n NUMBERS OF ELEMENTS ADJACENT TO ELEMENTS\n");
02081 for (i=0;i<ne;i++){
02082 fprintf (out,"\n element %5ld %ld",i,nadjelel[i]);
02083 }
02084
02085 fprintf (out,"\n\n\n LIST OF ELEMENTS ADJACENT TO ELEMENTS\n");
02086 for (i=0;i<ne;i++){
02087 fprintf (out,"\n element %4ld ",i);
02088 for (j=0;j<nadjelel[i];j++){
02089 fprintf (out," %ld",adjelel[i][j]);
02090 }
02091 }
02092 fprintf (out,"\n\n\n");
02093 }
02094
02095
02096 for (i=0;i<ne;i++){
02097 delete [] aux[i];
02098 }
02099 delete [] aux;
02100
02101 }
02102
02103
02104
02105
02106
02107
02108
02109
02110
02111 void gtopology::adjacnodes (FILE *out)
02112 {
02113 long i,j,k,l,m,min,prev,nne;
02114 long **auxadjnodnod;
02115 ivector nodes;
02116
02117
02118 if (nadjelnod==NULL)
02119 nadjelnod = new long [nn];
02120 memset (nadjelnod,0,nn*sizeof(long));
02121
02122
02123 for (i=0;i<ne;i++){
02124 nne = give_nne (i);
02125 reallocv (nne,nodes);
02126 give_nodes (i,nodes);
02127 for (j=0;j<nne;j++){
02128 nadjelnod[nodes[j]]++;
02129 }
02130 }
02131
02132
02133
02134
02135
02136
02137
02138
02139 if (adjelnod==NULL){
02140 adjelnod = new long* [nn];
02141 for (i=0;i<nn;i++){
02142 adjelnod[i] = new long [nadjelnod[i]];
02143 }
02144 }
02145 for (i=0;i<nn;i++){
02146 nadjelnod[i] = 0;
02147 }
02148
02149
02150 for (i=0;i<ne;i++){
02151 nne = give_nne (i);
02152 reallocv (nne,nodes);
02153 give_nodes (i,nodes);
02154 for (j=0;j<nne;j++){
02155 adjelnod[nodes[j]][nadjelnod[nodes[j]]++] = i;
02156 }
02157 }
02158
02159
02160
02161
02162
02163
02164
02165
02166
02167
02168
02169
02170 nadjnodnod = new long [nn];
02171 for (i=0;i<nn;i++){
02172 nadjnodnod[i]=0;
02173 for (j=0;j<nadjelnod[i];j++){
02174 k=adjelnod[i][j];
02175 nadjnodnod[i]+=give_nne (k);
02176 }
02177 }
02178
02179
02180 auxadjnodnod = new long* [nn];
02181 for (i=0;i<nn;i++){
02182 auxadjnodnod[i] = new long [nadjnodnod[i]];
02183 nadjnodnod[i]=0;
02184 }
02185
02186
02187 for (i=0;i<nn;i++){
02188 for (j=0;j<nadjelnod[i];j++){
02189 l=adjelnod[i][j];
02190 nne=give_nne (l);
02191 reallocv (nne,nodes);
02192 give_nodes (l,nodes);
02193 for (k=0;k<nne;k++){
02194 auxadjnodnod[i][nadjnodnod[i]]=nodes[k];
02195 nadjnodnod[i]++;
02196 }
02197 }
02198 }
02199
02200
02201
02202
02203
02204
02205
02206
02207
02208
02209
02210
02211
02212
02213
02214 for (i=0;i<nn;i++){
02215 prev=LONG_MAX;
02216 for (j=0;j<nadjnodnod[i];j++){
02217 min=LONG_MAX;
02218 for (k=j;k<nadjnodnod[i];k++){
02219 if (min>auxadjnodnod[i][k]){
02220 min=auxadjnodnod[i][k]; l=k;
02221 }
02222 }
02223 if (min==prev){
02224 nadjnodnod[i]--;
02225 auxadjnodnod[i][l]=auxadjnodnod[i][nadjnodnod[i]];
02226 j--;
02227 }
02228 else{
02229 m=auxadjnodnod[i][j];
02230 auxadjnodnod[i][j]=min;
02231 auxadjnodnod[i][l]=m;
02232 prev=min;
02233 }
02234 }
02235 }
02236
02237
02238
02239
02240
02241 adjnodnod = new long* [nn];
02242 for (i=0;i<nn;i++){
02243 adjnodnod[i] = new long [nadjnodnod[i]];
02244 for (j=0;j<nadjnodnod[i];j++){
02245 adjnodnod[i][j]=auxadjnodnod[i][j];
02246 }
02247 delete [] auxadjnodnod[i];
02248 }
02249 delete [] auxadjnodnod;
02250
02251
02252 if (out)
02253 {
02254 fprintf (out,"\n\n\n LIST OF NUMBERS OF NODES ADJACENT TO NODE\n");
02255 for (i=0;i<nn;i++)
02256 fprintf(out,"%ld %ld\n",i,nadjnodnod[i]);
02257 fprintf (out,"\n\n\n LIST OF NODES ADJACENT TO NODE\n");
02258 for (i=0;i<nn;i++){
02259 fprintf (out,"\n node %4ld ",i);
02260 for (j=0;j<nadjnodnod[i];j++){
02261 fprintf (out," %ld",adjnodnod[i][j]);
02262 }
02263 }
02264 fprintf (out,"\n\n\n");
02265 }
02266
02267 }
02268
02269
02270
02271
02272
02273
02274
02275
02276 void gtopology::adjacnodes_edge (FILE *out)
02277 {
02278 long i,j,k,l,m;
02279 long ned,nned,auxnedges;
02280 long *auxnned,*dupledges,*aux;
02281 long **auxedgenodes;
02282
02283
02284 auxnedges = 0;
02285 for(i = 0; i < ne; i++){
02286
02287 ned = give_ned(i);
02288 auxnedges+=ned;
02289 }
02290
02291
02292
02293 aux = new long[3];
02294
02295 auxnned = new long[auxnedges];
02296 auxedgenodes = new long*[auxnedges];
02297 dupledges = new long[auxnedges];
02298 k = 0;
02299 for(i = 0; i < ne; i++){
02300 ned = give_ned(i);
02301 nned = give_nned(i);
02302 for(j = 0; j < ned; j++ ){
02303 give_edge_nodes(i,j,aux);
02304 auxnned[k] = nned;
02305 auxedgenodes[k] = new long[nned];
02306 for(l = 0; l < nned; l++){
02307 auxedgenodes[k][l] = aux[l];
02308 }
02309
02310 if(nned == 2){
02311 if(auxedgenodes[k][0] > auxedgenodes[k][1]){
02312 l = auxedgenodes[k][1];
02313 auxedgenodes[k][1] = auxedgenodes[k][0];
02314 auxedgenodes[k][0] = l;
02315 }
02316 }
02317 if(nned == 3){
02318 if(auxedgenodes[k][0] > auxedgenodes[k][2]){
02319 l = auxedgenodes[k][2];
02320 auxedgenodes[k][2] = auxedgenodes[k][0];
02321 auxedgenodes[k][0] = l;
02322 }
02323 if(auxedgenodes[k][0] > auxedgenodes[k][1]){
02324 l = auxedgenodes[k][1];
02325 auxedgenodes[k][1] = auxedgenodes[k][0];
02326 auxedgenodes[k][0] = l;
02327 }
02328 if(auxedgenodes[k][1] > auxedgenodes[k][2]){
02329 l = auxedgenodes[k][2];
02330 auxedgenodes[k][2] = auxedgenodes[k][1];
02331 auxedgenodes[k][1] = l;
02332 }
02333 }
02334 dupledges[k] = 0;
02335 k++;
02336 }
02337 }
02338 delete []aux;
02339
02340
02341
02342
02343
02344
02345
02346
02347
02348
02349
02350
02351
02352
02353 for(i = 0; i < auxnedges; i++){
02354 for(j = i+1; j < auxnedges; j++){
02355 if(auxedgenodes[i][0] == auxedgenodes[j][0]){
02356 if(auxnned[i] == 2 && auxnned[j] == 2){
02357 if(auxedgenodes[i][1] == auxedgenodes[j][1]){
02358 dupledges[j]++;
02359
02360 }
02361 }
02362 if(auxnned[i] == 3 && auxnned[i] == 3){
02363 if(auxedgenodes[i][1] == auxedgenodes[j][1]){
02364 if(auxedgenodes[i][2] == auxedgenodes[j][2]){
02365 dupledges[j]++;
02366
02367 }
02368 }
02369 }
02370 }
02371 }
02372 }
02373
02374
02375
02376
02377
02378
02379
02380 nadjacnodesedge = new long[nn];
02381 for(i = 0; i < nn; i++){
02382 nadjacnodesedge[i] = 0;
02383 for(j = 0; j < auxnedges; j++){
02384 if(dupledges[j] == 0) {
02385 for(k = 0; k < auxnned[j]; k++){
02386 if(auxedgenodes[j][k] == i){
02387 auxedgenodes[j][k] = i;
02388 nadjacnodesedge[i]++;
02389 }
02390 }
02391 }
02392 }
02393 }
02394
02395
02396
02397
02398
02399
02400
02401
02402
02403
02404
02405
02406
02407 adjacnodesedge = new long*[nn];
02408 for(i = 0; i < nn; i++){
02409 adjacnodesedge[i] = new long[nadjacnodesedge[i]];
02410 m = 0;
02411 for(j = 0; j < auxnedges; j++){
02412 if(dupledges[j] == 0) {
02413 if(auxnned[j] == 2){
02414 for(k = 0; k < auxnned[j]; k++){
02415 if(auxedgenodes[j][k] == i){
02416 if(k == 0){
02417 adjacnodesedge[i][m] = auxedgenodes[j][1];
02418 m++;
02419 }
02420 if(k == 1){
02421 adjacnodesedge[i][m] = auxedgenodes[j][0];
02422 m++;
02423 }
02424 }
02425 }
02426 }
02427 if(auxnned[j] == 3){
02428 for(k = 0; k < auxnned[j]; k++){
02429 if(auxedgenodes[j][k] == i){
02430 if(k == 0){
02431 adjacnodesedge[i][m] = auxedgenodes[j][1];
02432 m++;
02433 }
02434 if(k == 1){
02435 adjacnodesedge[i][m] = auxedgenodes[j][0];
02436 m++;
02437 adjacnodesedge[i][m] = auxedgenodes[j][2];
02438 m++;
02439 }
02440 if(k == 2){
02441 adjacnodesedge[i][m] = auxedgenodes[j][1];
02442 m++;
02443 }
02444 }
02445 }
02446 }
02447 }
02448 }
02449 }
02450
02451
02452
02453
02454
02455
02456
02457
02458
02459
02460
02461 for(i = 0; i < nn; i++){
02462 if(nadjacnodesedge[i] == 2){
02463 if(adjacnodesedge[i][0] > adjacnodesedge[i][1]){
02464 j = adjacnodesedge[i][1];
02465 adjacnodesedge[i][1] = adjacnodesedge[i][0];
02466 adjacnodesedge[i][0] = j;
02467 }
02468 }
02469 else{
02470 for(j = nadjacnodesedge[i] - 1; j > 0; j--){
02471 for(k = 0; k < j; k++){
02472 if(adjacnodesedge[i][k] > adjacnodesedge[i][k+1]){
02473 m = adjacnodesedge[i][k];
02474 adjacnodesedge[i][k] = adjacnodesedge[i][k+1];
02475 adjacnodesedge[i][k+1] = m;
02476 }
02477 }
02478 }
02479 }
02480 }
02481
02482
02483
02484 for(i = 0; i < nn; i++){
02485 fprintf(out,"node %ld has %ld neighbours: ",i,nadjacnodesedge[i]);
02486 for(j = 0; j < nadjacnodesedge[i]; j++){
02487 fprintf(out," %ld",adjacnodesedge[i][j]);
02488 }
02489 fprintf(out,"\n");
02490 }
02491
02492
02493 for(i = 0; i < auxnedges; i++){
02494 delete []auxedgenodes[i];
02495 }
02496 delete []auxnned;
02497 delete []auxedgenodes;
02498 delete []dupledges;
02499
02500 }
02501
02502
02503
02504
02505
02506
02507
02508
02509
02510
02511
02512 void gtopology::backup_codnum (void)
02513 {
02514 long i, j, ndofn;
02515
02516 for (i = 0; i < nn; i++)
02517 {
02518 ndofn = give_ndofn (i);
02519 if (ndofn < 0)
02520 ndofn = -ndofn;
02521
02522 bckcn[i] = new long[ndofn];
02523 for (j = 0; j < ndofn; j++)
02524 bckcn[i][j]=give_dof (i,j);
02525 }
02526 }
02527
02528
02529
02530
02531
02532
02533
02534
02535 void gtopology::restore_codnum (void)
02536 {
02537 long i, j, ndofn;
02538
02539 for (i = 0; i < nn; i++)
02540 {
02541 ndofn = give_ndofn (i);
02542 if (ndofn < 0)
02543 ndofn = -ndofn;
02544 for (j = 0; j < ndofn; j++)
02545 gnodes[i].cn[j] = bckcn[i][j];
02546 }
02547 }
02548
02549 void gtopology::write_nodes(FILE *out)
02550 {
02551 fprintf(out, "NODES\n");
02552 for (long i = 0; i < nn; i++)
02553 fprintf(out, "%ld %e %e %e\n", i, gnodes[i].x, gnodes[i].y, gnodes[i].z);
02554 }
02555
02556
02557
02558
02559
02560
02561
02562
02563 void gtopology::unodelnode ()
02564 {
02565 long i,j;
02566 unln = new long [nn];
02567 unnl = new long [nn];
02568 for (i=0;i<nln;i++){
02569 for (j=0;j<lgnodes[i].nl;j++){
02570 unln[lgnodes[i].nodes[j]]=i;
02571 unnl[lgnodes[i].nodes[j]]=j;
02572 }
02573 }
02574 }
02575
02576
02577
02578
02579
02580
02581 void gtopology::initiate_elemmult ()
02582 {
02583 long i,nid,lnid;
02584
02585 for (i=0;i<ne;i++){
02586 nid=gelements[i].nodes[0];
02587 lnid=unln[nid];
02588 gelements[i].nmult=give_nmult (lnid);
02589 }
02590 }
02591
02592
02593
02594
02595
02596
02597 void gtopology::comp_domain_sizes ()
02598 {
02599 long i;
02600 double x,y,z,minx,miny,minz,maxx,maxy,maxz;
02601
02602 domsizes = new double [3];
02603
02604 minx=gnodes[0].x; maxx=gnodes[0].x;
02605 miny=gnodes[0].y; maxy=gnodes[0].y;
02606 minz=gnodes[0].z; maxz=gnodes[0].z;
02607
02608 for (i=1;i<nn;i++){
02609 x=gnodes[i].x;
02610 y=gnodes[i].y;
02611 z=gnodes[i].z;
02612
02613 if (minx>x) minx=x;
02614 if (maxx<x) maxx=x;
02615
02616 if (miny>y) miny=y;
02617 if (maxy<y) maxy=y;
02618
02619 if (minz>z) minx=z;
02620 if (maxz<z) maxx=z;
02621 }
02622
02623 domsizes[0]=maxx-minx;
02624 domsizes[1]=maxy-miny;
02625 domsizes[2]=maxz-minz;
02626 }
02627
02628
02629
02630
02631
02632
02633 void gtopology::give_domain_sizes (double *sizes)
02634 {
02635 sizes[0]=domsizes[0];
02636 sizes[1]=domsizes[1];
02637 sizes[2]=domsizes[2];
02638 }
02639
02640
02641
02642
02643
02644
02645
02646
02647
02648
02649
02650 void gtopology::read_gf (XFILE *in)
02651 {
02652 long i,j;
02653 char emsg[200];
02654 xfscanf (in,"%k%ld","num_gfunct", &ngf);
02655 gf = new gfunct [ngf];
02656
02657 for (i=0;i<ngf;i++){
02658 xfscanf (in,"%k%ld", "gf_id",&j);
02659 if ((j < 1) || (j > ngf))
02660 {
02661 sprintf(emsg, "general function index %ld out of range <1,%ld>\n"
02662 " at line=%ld, col=%ld, file %s", j, ngf, in->line, in->col, in->fname);
02663 print_err(emsg, __FILE__, __LINE__, __func__);
02664 abort();
02665 }
02666 j--;
02667 gf[j].read (in);
02668 }
02669 }
02670
02671
02672
02673
02674
02675
02676
02677
02678
02679 void gtopology::print_gf (FILE *out)
02680 {
02681 long i;
02682 fprintf (out,"\n\n%ld\n",ngf);
02683
02684 for (i=0;i<ngf;i++){
02685 fprintf (out,"%ld\n",i+1);
02686 gf[i].print (out);
02687 }
02688 }
02689
02690
02691
02692
02693
02694
02695 void gtopology::lneso_init ()
02696 {
02697 long i;
02698
02699
02700 leso = new long [ne];
02701
02702 for (i=0;i<ne;i++){
02703 leso[i]=1;
02704 }
02705
02706
02707 lnso = new long [nn];
02708
02709 for (i=0;i<nn;i++){
02710 lnso[i]=1;
02711 }
02712 }
02713
02714
02715
02716
02717
02718
02719
02720 void gtopology::auxinf_init ()
02721 {
02722 long i;
02723
02724 for (i=0;i<ne;i++){
02725 gelements[i].auxinf=0;
02726 }
02727 }
02728
02729
02730
02731
02732
02733
02734
02735
02736
02737
02738
02739 void gtopology::update_auxinf()
02740 {
02741 long i;
02742
02743 for (i=0;i<ne;i++)
02744 gelements[i].auxinf = leso[i];
02745 }
02746
02747
02748
02749
02750
02751
02752
02753
02754
02755
02756
02757
02758
02759
02760
02761
02762 void gtopology::update_elem (double time)
02763 {
02764 long i,j,k;
02765
02766 for (i=0;i<ne;i++){
02767
02768 j=gelements[i].tgf;
02769
02770 k=gf[j].getval_long (time);
02771
02772 leso[i]=k;
02773 }
02774 }
02775
02776
02777
02778
02779
02780
02781
02782
02783
02784
02785
02786 void gtopology::update_nodes ()
02787 {
02788 long i,j,nne, nmne;
02789 ivector nod;
02790
02791 for (i=0;i<nn;i++){
02792 lnso[i]=0;
02793 }
02794
02795 for (i=0;i<ne;i++){
02796 if (leso[i]==1){
02797
02798
02799
02800 nne = give_nne (i);
02801 reallocv (nne,nod);
02802
02803 give_nodes (i,nod);
02804 for (j=0;j<nne;j++){
02805 lnso[nod[j]]=1;
02806 }
02807 nmne = give_nmne(i);
02808 if (nmne > 0)
02809 {
02810
02811 nne = give_original_nne(i);
02812 reallocv(nne, nod);
02813 give_original_nodes(i,nod);
02814 for (j=0;j<nne;j++)
02815 {
02816 if (gnodes[nod[j]].mnodes)
02817 lnso[nod[j]]=1;
02818 }
02819 }
02820 }
02821 }
02822 }
02823
02824
02825
02826
02827
02828
02829
02830
02831
02832
02833
02834
02835
02836
02837
02838 void gtopology::remove_nodes(long *sn)
02839 {
02840 long i;
02841 ivector nod;
02842
02843 for (i=0;i<nn;i++)
02844 {
02845 if (sn[i])
02846 lnso[i]=0;
02847 }
02848 }
02849
02850
02851
02852
02853
02854
02855
02856
02857
02858
02859
02860
02861
02862
02863
02864
02865
02866 long gtopology::search_changed_elem(double time, long &nae)
02867 {
02868 long i, j, k, l;
02869 long nce;
02870
02871
02872 nae=0;
02873
02874 nce = 0;
02875
02876 for (i=0;i<ne;i++)
02877 {
02878
02879 j=gelements[i].tgf;
02880
02881 k=gf[j].getval_long(time);
02882
02883 l=gelements[i].auxinf;
02884
02885 leso[i] = k;
02886
02887 if (k==l)
02888 continue;
02889 else
02890 {
02891 if (k)
02892 nae++;
02893 nce++;
02894 }
02895 }
02896 return nce;
02897 }
02898
02899
02900
02901
02902
02903
02904
02905
02906
02907
02908
02909
02910
02911
02912
02913
02914 long gtopology::search_iface_nodes(long *ifn)
02915 {
02916 long i, j, nifn, nne;
02917 ivector nod;
02918
02919 memset(ifn, 0, sizeof(*ifn)*nn);
02920
02921
02922 for (i=0;i<ne;i++)
02923 {
02924 if ((leso[i]) && (gelements[i].auxinf))
02925 {
02926 nne = give_nne(i);
02927 reallocv(nne,nod);
02928
02929 give_nodes(i,nod);
02930
02931 for(j=0; j<nne; j++)
02932 {
02933 if (ifn[nod[j]] == 0)
02934 ifn[nod[j]] = 1;
02935 }
02936 }
02937 }
02938
02939 for (i=0;i<ne;i++)
02940 {
02941 if ((leso[i]) && (gelements[i].auxinf == 0))
02942 {
02943 nne = give_nne(i);
02944 reallocv(nne,nod);
02945
02946 give_nodes(i,nod);
02947 for(j=0; j<nne; j++)
02948 {
02949 if (ifn[nod[j]] == 1)
02950 ifn[nod[j]] = 2;
02951 }
02952 }
02953 }
02954 nifn = 0;
02955 for (i=0; i<nn; i++)
02956 {
02957 if (ifn[i] < 2)
02958 ifn[i] = 0;
02959 else
02960 {
02961 ifn[i] = 1;
02962 nifn++;
02963 }
02964 }
02965 return nifn;
02966 }
02967
02968
02969
02970
02971
02972
02973
02974
02975
02976
02977
02978
02979
02980
02981 long gtopology::search_newelem (double time,double prev_time)
02982 {
02983 long i,j,k,l,nae;
02984
02985
02986 nae=0;
02987
02988 if (time>prev_time){
02989 for (i=0;i<ne;i++){
02990
02991 j=gelements[i].tgf;
02992
02993 k=gf[j].getval_long (time);
02994
02995 l=gf[j].getval_long (prev_time);
02996
02997 if (k==l){
02998 leso[i]=0;
02999 }
03000 else{
03001 leso[i]=1;
03002 nae++;
03003 }
03004 }
03005 }
03006 else{
03007
03008
03009 for (i=0;i<ne;i++){
03010
03011 j=gelements[i].tgf;
03012
03013 k=gf[j].getval_long (time);
03014
03015 if (k==0){
03016 leso[i]=0;
03017 }
03018 else{
03019 leso[i]=1;
03020 nae++;
03021 }
03022 }
03023 }
03024
03025 return nae;
03026 }
03027
03028
03029
03030
03031
03032
03033
03034
03035
03036
03037
03038
03039 long gtopology::search_newdofs (double time,double prev_time)
03040 {
03041 long i,nch;
03042 long *plnso = new long[nn];
03043
03044 memcpy(plnso, lnso, sizeof(*plnso)*nn);
03045
03046
03047
03048 update_nodes();
03049
03050
03051 nch=0;
03052 for (i=0;i<nn;i++){
03053 nch += gnodes[i].search_changed_dofs (gf,time,prev_time,lnso[i],plnso[i]);
03054 }
03055
03056
03057 if (nch > 0)
03058 cnstate=0;
03059
03060 delete [] plnso;
03061 return nch;
03062 }
03063
03064
03065
03066
03067
03068
03069
03070
03071
03072
03073
03074
03075
03076 long gtopology::search_changed_dofs (double time,double prev_time)
03077 {
03078 long i,nch;
03079 long *plnso = new long[nn];
03080
03081 memcpy(plnso, lnso, sizeof(*plnso)*nn);
03082
03083
03084
03085 update_nodes();
03086
03087
03088 nch=0;
03089 for (i=0;i<nn;i++){
03090 nch += gnodes[i].search_changed_dofs(gf,time,prev_time,lnso[i],plnso[i]);
03091 }
03092
03093
03094 if (nch > 0)
03095 cnstate=0;
03096
03097 delete [] plnso;
03098 return nch;
03099 }
03100
03101
03102
03103
03104
03105
03106
03107
03108
03109
03110
03111 void gtopology::switch_new_elem()
03112 {
03113 long i;
03114 for (i=0; i<ne; i++)
03115 {
03116 if (leso[i] && (gelements[i].auxinf == 0))
03117 continue;
03118 else
03119 leso[i] = 0;
03120 }
03121 }
03122
03123
03124
03125
03126
03127
03128
03129
03130
03131
03132
03133
03134
03135 void gtopology::update_active_dofs (double time)
03136 {
03137 long i;
03138
03139 for (i=0;i<nn;i++){
03140 if (lnso[i])
03141 gnodes[i].update_dofs (gf,time,lnso[i]);
03142 else
03143 {
03144 if (gnodes[i].ndofn > 0)
03145 gnodes[i].clear_dof();
03146 }
03147 }
03148
03149 for (i=0;i<ne;i++)
03150 {
03151 if (gelements[i].master_nodes!=NULL)
03152 {
03153
03154 gelements[i].cne=0;
03155 }
03156 }
03157
03158
03159 cnstate=0;
03160 }
03161
03162
03163
03164
03165
03166
03167
03168
03169
03170
03171
03172
03173 void gtopology::clear_intf_dofs(long *ifn)
03174 {
03175 long i;
03176
03177 for (i=0;i<nn;i++)
03178 {
03179 if (ifn[i])
03180 gnodes[i].clear_dof();
03181 }
03182
03183
03184 cnstate=0;
03185 }
03186
03187
03188
03189
03190
03191
03192
03193
03194
03195
03196
03197
03198 void gtopology::update_dofs (double time)
03199 {
03200 long i;
03201
03202 for (i=0;i<nn;i++)
03203 gnodes[i].update_dofs (gf,time,lnso[i]);
03204
03205
03206 cnstate=0;
03207 }
03208
03209
03210
03211
03212
03213
03214
03215
03216
03217
03218
03219
03220
03221
03222
03223
03224
03225
03226
03227
03228
03229
03230
03231 void gtopology::print_time_functions (FILE *kon,double time)
03232 {
03233 long i;
03234
03235
03236
03237
03238
03239
03240 fprintf (kon,"\n\n vypis v case %le\n",time);
03241 for (i=0;i<ngf;i++){
03242
03243 fprintf (kon,"%ld %ld\n",i,gf[i].getval_long (time));
03244 }
03245
03246
03247 }
03248
03249
03250
03251
03252
03253
03254
03255
03256
03257
03258
03259
03260
03261
03262
03263
03264 void gtopology::end_nodes (FILE *out)
03265 {
03266 long i,j,k,nsn,nene;
03267 long *nodes,*sn;
03268 long *av;
03269
03270 nodes = new long[2];
03271
03272
03273 av = new long [nn];
03274 for (i=0;i<nn;i++){
03275 av[i]=0;
03276 }
03277
03278
03279 for (i=0;i<ne;i++){
03280
03281 nene=give_nen (i);
03282
03283 if (nene>0){
03284
03285 give_end_nodes (i,nodes);
03286
03287
03288 for (j=0;j<2;j++){
03289
03290
03291
03292
03293 av[nodes[j]]++;
03294 }
03295 }
03296 }
03297
03298
03299 fprintf (out,"\n\n\n kontrola koncovych bodu \n\n");
03300 for (i=0;i<nn;i++){
03301 fprintf (out,"\n node %4ld av %4ld",i,av[i]);
03302 }
03303 fprintf (out,"\n\n\n");
03304
03305
03306 nsn=0;
03307 for (i=0;i<nn;i++){
03308 if (av[i]==1)
03309 nsn++;
03310 }
03311
03312
03313 sn = new long [nsn];
03314 nsn=0;
03315 for (i=0;i<nn;i++){
03316 if (av[i]==1){
03317 sn[nsn]=i;
03318 nsn++;
03319 }
03320 }
03321
03322
03323 nen=0;
03324 double threshold=1.0e-7;
03325 for (i=0;i<nsn;i++){
03326 for (j=i+1;j<nsn;j++){
03327 if (fabs(gnodes[sn[j]].x-gnodes[sn[i]].x)<threshold){
03328 if (fabs(gnodes[sn[j]].y-gnodes[sn[i]].y)<threshold){
03329 if (fabs(gnodes[sn[j]].z-gnodes[sn[i]].z)<threshold){
03330 nen++;
03331 }
03332 }
03333 }
03334 }
03335 }
03336
03337 fprintf (stdout,"\n pocet zajimavych koncovych bodu %ld",nen);
03338 fprintf (out,"\n\n\n pocet zajimavych koncovych bodu %ld \n\n",nen);
03339
03340 endnodes = new endnode [nen];
03341 k=0;
03342 for (i=0;i<nsn;i++){
03343 for (j=i+1;j<nsn;j++){
03344 if (fabs(gnodes[sn[j]].x-gnodes[sn[i]].x)<threshold){
03345 if (fabs(gnodes[sn[j]].y-gnodes[sn[i]].y)<threshold){
03346 if (fabs(gnodes[sn[j]].z-gnodes[sn[i]].z)<threshold){
03347 endnodes[k].nn=2;
03348 endnodes[k].nm=2;
03349 endnodes[k].fn=sn[i];
03350 endnodes[k].ln=sn[j];
03351 k++;
03352 }
03353 }
03354 }
03355 }
03356 }
03357
03358 delete [] nodes;
03359 delete [] sn;
03360 delete [] av;
03361 }
03362
03363
03364
03365
03366
03367
03368
03369
03370 void gtopology::endnodes_auxprint (FILE *out)
03371 {
03372 long i;
03373
03374 fprintf (out,"\n\n number of end nodes %ld",nen);
03375
03376 for (i=0;i<nen;i++){
03377 fprintf (out,"\n\n end node number %6ld",i);
03378 fprintf (out,"\n number of nodes %ld",endnodes[i].nn);
03379 fprintf (out,"\n number of Lagrange multipliers %ld",endnodes[i].nm);
03380 fprintf (out,"\n node number 1 %ld",endnodes[i].fn);
03381 fprintf (out,"\n node number 2 %ld",endnodes[i].ln);
03382 }
03383 }
03384
03385
03386
03387
03388
03389
03390
03391 void gtopology::alloc_endnode_cn ()
03392 {
03393 long i,nid,ndofn;
03394
03395 for (i=0;i<nen;i++){
03396
03397 nid=endnodes[i].fn;
03398
03399 ndofn=give_ndofn (nid);
03400
03401
03402 endnodes[i].alloc_cnm (ndofn);
03403 }
03404
03405 }
03406
03407
03408
03409
03410
03411
03412
03413
03414 void gtopology::endnodes_localization (gmatrix *gm)
03415 {
03416 long i,nm,*ncn1,*ncn2,*mcn;
03417
03418 for (i=0;i<nen;i++){
03419
03420 nm=endnodes[i].ndofn;
03421 ncn1 = new long [nm];
03422 ncn2 = new long [nm];
03423 mcn = new long [nm];
03424
03425
03426 give_endnode_code_numbers (i,ncn1,ncn2,mcn);
03427
03428
03429 gm->mult_localize (nm,ncn1,ncn2,mcn);
03430
03431 delete [] mcn;
03432 delete [] ncn2;
03433 delete [] ncn1;
03434 }
03435 }
03436
03437
03438
03439
03440
03441
03442 void gtopology::enodes_allelem (FILE *out)
03443 {
03444 long i,k,l,fn,nne,stop;
03445 ivector nod;
03446
03447
03448 for (i=0;i<nen;i++){
03449
03450 endnodes[i].adjel = new long [2];
03451
03452
03453 fn=endnodes[i].fn;
03454
03455
03456
03457 stop=0;
03458 for (k=0;k<ne;k++){
03459
03460 nne = give_original_nne (k);
03461
03462 allocv (nne,nod);
03463
03464 give_original_nodes (k,nod);
03465
03466
03467 for (l=0;l<nne;l++){
03468 if (nod[l]==fn){
03469 endnodes[i].adjel[0]=k;
03470 stop=1;
03471 break;
03472 }
03473 }
03474 destrv (nod);
03475 if (stop==1)
03476 break;
03477 }
03478
03479
03480
03481 fn=endnodes[i].ln;
03482
03483
03484
03485 stop=0;
03486 for (k=0;k<ne;k++){
03487
03488 nne = give_nne (k);
03489
03490 allocv (nne,nod);
03491
03492 give_original_nodes (k,nod);
03493
03494
03495 for (l=0;l<nne;l++){
03496 if (nod[l]==fn){
03497 endnodes[i].adjel[1]=k;
03498 stop=1;
03499 break;
03500 }
03501 }
03502 destrv (nod);
03503 if (stop==1)
03504 break;
03505 }
03506 }
03507 }
03508
03509
03510
03511
03512
03513
03514
03515
03516
03517
03518
03519
03520
03521
03522
03523
03524
03525
03526
03527
03528
03529
03530
03531
03532
03533
03534
03535
03536
03537
03538
03539
03540
03541
03542
03543
03544
03545
03546
03547
03548
03549
03550
03551
03552
03553
03554
03555
03556
03557
03558
03559
03560
03561
03562 void gtopology::edges_mat (long matrad,FILE *out)
03563 {
03564 long i,j,k,l,m,elid,nsed,ned,nedk,nned,nnedk,n1,n2,n3,n4,aux,mat1,mat2,stopp;
03565 long *nodes,*nodesk;
03566
03567
03568 nsed=0;
03569
03570
03571 for (i=0;i<ne;i++){
03572
03573 mat1=gelements[i].auxinf;
03574
03575 ned=give_ned (i);
03576
03577 nned=give_nned (i);
03578
03579 nodes=new long [nned];
03580
03581
03582 for (j=0;j<ned;j++){
03583
03584 give_edge_nodes (i,j,nodes);
03585
03586 n1=nodes[0];
03587 n2=nodes[nned-1];
03588 if (n1>n2){
03589 aux=n2;
03590 n2=n1;
03591 n1=aux;
03592 }
03593
03594 stopp=0;
03595
03596
03597 for (k=0;k<nadjelel[i];k++){
03598
03599 elid=adjelel[i][k];
03600
03601
03602 mat2=gelements[elid].auxinf;
03603
03604 nedk=give_ned (elid);
03605
03606 nnedk=give_nned (elid);
03607
03608 nodesk=new long [nnedk];
03609
03610
03611 for (l=0;l<nedk;l++){
03612
03613 give_edge_nodes (elid,l,nodesk);
03614
03615 n3=nodesk[0];
03616 n4=nodesk[nnedk-1];
03617 if (n3>n4){
03618 aux=n4;
03619 n4=n3;
03620 n3=aux;
03621 }
03622
03623 if (n1==n3 && n2==n4){
03624 if ((mat1==matrad && mat2!=matrad) || (mat1!=matrad && mat2==matrad)){
03625 nsed++;
03626
03627
03628 stopp=1;
03629 }
03630 }
03631 if (stopp==1){
03632 break;
03633 }
03634 }
03635 delete [] nodesk;
03636 if (stopp==1){
03637 break;
03638 }
03639 }
03640 }
03641 delete [] nodes;
03642 }
03643
03644
03645
03646 nged=nsed/2;
03647
03648
03649 gedges = new gedge [nged];
03650
03651
03652 nsed=0;
03653
03654 for (i=0;i<ne;i++){
03655
03656 mat1=gelements[i].auxinf;
03657
03658 ned=give_ned (i);
03659
03660 nned=give_nned (i);
03661
03662 nodes=new long [nned];
03663
03664
03665 for (j=0;j<ned;j++){
03666
03667 give_edge_nodes (i,j,nodes);
03668
03669 n1=nodes[0];
03670 n2=nodes[nned-1];
03671 if (n1>n2){
03672
03673
03674 aux=n2;
03675 n2=n1;
03676 n1=aux;
03677 }
03678
03679 stopp=0;
03680
03681
03682 for (k=0;k<nadjelel[i];k++){
03683
03684 elid=adjelel[i][k];
03685
03686
03687 mat2=gelements[elid].auxinf;
03688
03689 nedk=give_ned (elid);
03690
03691 nnedk=give_nned (elid);
03692
03693 nodesk=new long [nnedk];
03694
03695
03696 for (l=0;l<nedk;l++){
03697
03698 give_edge_nodes (elid,l,nodesk);
03699
03700 n3=nodesk[0];
03701 n4=nodesk[nnedk-1];
03702 if (n3>n4){
03703
03704
03705 aux=n4;
03706 n4=n3;
03707 n3=aux;
03708 }
03709
03710 if (n1==n3 && n2==n4){
03711
03712 if ((mat1==matrad && mat2!=matrad) || (mat1!=matrad && mat2==matrad)){
03713
03714
03715 for (m=0;m<nsed;m++){
03716 if (gedges[m].fn==n1 && gedges[m].ln==n2){
03717
03718 stopp=1;
03719 break;
03720 }
03721 }
03722 if (stopp==0){
03723
03724
03725 gedges[nsed].fn=n1;
03726 gedges[nsed].ln=n2;
03727 gedges[nsed].nn=4;
03728 gedges[nsed].nm=2;
03729
03730 gedges[nsed].nlist = new long [4];
03731 gedges[nsed].nlist[0]=n1;
03732 gedges[nsed].nlist[1]=n1;
03733 gedges[nsed].nlist[2]=n2;
03734 gedges[nsed].nlist[3]=n2;
03735
03736 gedges[nsed].adjel = new long [2];
03737 gedges[nsed].adjel[0]=i;
03738 gedges[nsed].adjel[1]=elid;
03739
03740 nsed++;
03741 stopp=1;
03742 }
03743 }
03744 }
03745 if (stopp==1){
03746 break;
03747 }
03748 }
03749 delete [] nodesk;
03750 if (stopp==1){
03751 break;
03752 }
03753 }
03754 }
03755 delete [] nodes;
03756 }
03757
03758 }
03759
03760
03761
03762
03763
03764
03765
03766
03767
03768
03769
03770
03771
03772
03773 void gtopology::edges (FILE *out)
03774 {
03775 long i,j,k,ned,nned,nsed;
03776 long *nodes,*fn,*ln;
03777 long **av;
03778
03779
03780
03781 av = new long* [nn];
03782 for (i=0;i<nn;i++){
03783 av[i] = new long [nadjnodnod[i]];
03784 for (j=0;j<nadjnodnod[i];j++){
03785 av[i][j]=0;
03786 }
03787 }
03788
03789
03790 for (i=0;i<ne;i++){
03791
03792 ned=give_ned (i);
03793
03794 nned=give_nned (i);
03795
03796 nodes=new long [nned];
03797
03798
03799 for (j=0;j<ned;j++){
03800
03801 give_edge_nodes (i,j,nodes);
03802
03803
03804
03805
03806 for (k=0;k<nadjnodnod[nodes[0]];k++){
03807 if (adjnodnod[nodes[0]][k]==nodes[nned-1]){
03808 av[nodes[0]][k]++;
03809 break;
03810 }
03811 }
03812 for (k=0;k<nadjnodnod[nodes[nned-1]];k++){
03813 if (adjnodnod[nodes[nned-1]][k]==nodes[0]){
03814 av[nodes[nned-1]][k]++;
03815 break;
03816 }
03817 }
03818 }
03819
03820 delete [] nodes;
03821 }
03822
03823
03824 fprintf (out,"\n\n\n kontrola hran \n\n");
03825 for (i=0;i<nn;i++){
03826 fprintf (out,"\n node %4ld ",i);
03827 for (j=0;j<nadjnodnod[i];j++){
03828 fprintf (out," %4ld : %2ld",adjnodnod[i][j],av[i][j]);
03829 }
03830 }
03831 fprintf (out,"\n\n\n");
03832
03833
03834
03835 nsed=0;
03836 for (i=0;i<nn;i++){
03837 for (j=0;j<nadjnodnod[i];j++){
03838 if (av[i][j]==1)
03839 nsed++;
03840 }
03841 }
03842 nsed/=2;
03843
03844
03845 fn = new long [nsed];
03846
03847 ln = new long [nsed];
03848
03849 k=0;
03850 for (i=0;i<nn;i++){
03851 for (j=0;j<nadjnodnod[i];j++){
03852 if (av[i][j]==1){
03853 if (i<adjnodnod[i][j]){
03854 fn[k]=i;
03855 ln[k]=adjnodnod[i][j];
03856 k++;
03857 }
03858 }
03859 }
03860 }
03861
03862 fprintf (out,"\n\n\n druha kontrola hran \n\n");
03863 for (i=0;i<nsed;i++){
03864 fprintf (out,"\n fn %5ld ln %5ld",fn[i]+1,ln[i]+1);
03865 }
03866
03867
03868
03869
03870
03871 nged=0;
03872 double threshold=1.0e-7;
03873 for (i=0;i<nsed;i++){
03874 for (j=i+1;j<nsed;j++){
03875 if (fabs(gnodes[fn[j]].x-gnodes[fn[i]].x)<threshold && fabs(gnodes[ln[j]].x-gnodes[ln[i]].x)<threshold){
03876 if (fabs(gnodes[fn[j]].y-gnodes[fn[i]].y)<threshold && fabs(gnodes[ln[j]].y-gnodes[ln[i]].y)<threshold){
03877 if (fabs(gnodes[fn[j]].z-gnodes[fn[i]].z)<threshold && fabs(gnodes[ln[j]].z-gnodes[ln[i]].z)<threshold){
03878 nged++;
03879 }
03880 }
03881 }
03882 if (fabs(gnodes[fn[j]].x-gnodes[ln[i]].x)<threshold && fabs(gnodes[ln[j]].x-gnodes[fn[i]].x)<threshold){
03883 if (fabs(gnodes[fn[j]].y-gnodes[ln[i]].y)<threshold && fabs(gnodes[ln[j]].y-gnodes[fn[i]].y)<threshold){
03884 if (fabs(gnodes[fn[j]].z-gnodes[ln[i]].z)<threshold && fabs(gnodes[ln[j]].z-gnodes[fn[i]].z)<threshold){
03885 nged++;
03886 }
03887 }
03888 }
03889 }
03890 }
03891
03892 fprintf (out,"\n\n\n pocet zajimavych hran %ld \n\n",nged);
03893
03894 gedges = new gedge [nged];
03895 k=0;
03896 for (i=0;i<nsed;i++){
03897 for (j=i+1;j<nsed;j++){
03898 if (fabs(gnodes[fn[j]].x-gnodes[fn[i]].x)<threshold && fabs(gnodes[ln[j]].x-gnodes[ln[i]].x)<threshold){
03899 if (fabs(gnodes[fn[j]].y-gnodes[fn[i]].y)<threshold && fabs(gnodes[ln[j]].y-gnodes[ln[i]].y)<threshold){
03900 if (fabs(gnodes[fn[j]].z-gnodes[fn[i]].z)<threshold && fabs(gnodes[ln[j]].z-gnodes[ln[i]].z)<threshold){
03901 gedges[k].nn=4;
03902 gedges[k].nm=2;
03903 gedges[k].fn=fn[i];
03904 gedges[k].ln=ln[i];
03905 gedges[k].nlist = new long [4];
03906 gedges[k].nlist[0]=fn[i];
03907 gedges[k].nlist[1]=fn[j];
03908 gedges[k].nlist[2]=ln[i];
03909 gedges[k].nlist[3]=ln[j];
03910 k++;
03911 }
03912 }
03913 }
03914 if (fabs(gnodes[fn[j]].x-gnodes[ln[i]].x)<threshold && fabs(gnodes[ln[j]].x-gnodes[fn[i]].x)<threshold){
03915 if (fabs(gnodes[fn[j]].y-gnodes[ln[i]].y)<threshold && fabs(gnodes[ln[j]].y-gnodes[fn[i]].y)<threshold){
03916 if (fabs(gnodes[fn[j]].z-gnodes[ln[i]].z)<threshold && fabs(gnodes[ln[j]].z-gnodes[fn[i]].z)<threshold){
03917 gedges[k].nn=4;
03918 gedges[k].nm=2;
03919 gedges[k].fn=fn[i];
03920 gedges[k].ln=ln[i];
03921 gedges[k].nlist = new long [4];
03922 gedges[k].nlist[0]=fn[i];
03923 gedges[k].nlist[1]=ln[j];
03924 gedges[k].nlist[2]=ln[i];
03925 gedges[k].nlist[3]=fn[j];
03926 k++;
03927 }
03928 }
03929 }
03930 }
03931 }
03932
03933
03934
03935
03936
03937 delete [] ln;
03938 delete [] fn;
03939 for (i=0;i<nn;i++){
03940 delete [] av[i];
03941 }
03942 delete [] av;
03943
03944 }
03945
03946
03947
03948
03949
03950
03951 void gtopology::edge_dirvect ()
03952 {
03953 long i;
03954
03955 for (i=0;i<nged;i++){
03956 gedges[i].direction_vector (gnodes);
03957 }
03958 }
03959
03960
03961
03962
03963
03964
03965 void gtopology::edge_normvect ()
03966 {
03967 long i;
03968
03969 for (i=0;i<nged;i++){
03970 gedges[i].normal_vector (gnodes);
03971 }
03972 }
03973
03974
03975
03976
03977
03978
03979 void gtopology::edgenode_sorting ()
03980 {
03981 long i,j,nc;
03982
03983 nc=0;
03984 for (i=0;i<nged;i++){
03985 if (gedges[i].nm != 2 || gedges[i].nn != 4){
03986 print_err(" wrong number of node multiplicity on edge", __FILE__, __LINE__, __func__);
03987 abort ();
03988 }
03989 if (gedges[i].nlist[0]>gedges[i].nlist[1]){
03990 j=gedges[i].nlist[1];
03991 gedges[i].nlist[1]=gedges[i].nlist[0];
03992 gedges[i].nlist[0]=j;
03993 nc++;
03994 }
03995 if (gedges[i].nlist[2]>gedges[i].nlist[3]){
03996 j=gedges[i].nlist[3];
03997 gedges[i].nlist[3]=gedges[i].nlist[2];
03998 gedges[i].nlist[2]=j;
03999 nc++;
04000 }
04001 }
04002
04003 fprintf (stdout,"\n number of corrections in function edgenode_sorting is %ld",nc);
04004 }
04005
04006
04007
04008
04009
04010
04011
04012 void gtopology::prev_next_edges ()
04013 {
04014 long i,j;
04015 long ifn1,ifn2,iln1,iln2,jfn1,jfn2,jln1,jln2;
04016
04017 for (i=0;i<nged;i++){
04018 if (gedges[i].nm != 2 || gedges[i].nn != 4){
04019 print_err(" wrong number of node multiplicity on edge", __FILE__, __LINE__, __func__);
04020 abort ();
04021 }
04022
04023 ifn1=gedges[i].nlist[0];
04024 ifn2=gedges[i].nlist[1];
04025 iln1=gedges[i].nlist[2];
04026 iln2=gedges[i].nlist[3];
04027
04028 for (j=i+1;j<nged;j++){
04029 jfn1=gedges[j].nlist[0];
04030 jfn2=gedges[j].nlist[1];
04031 jln1=gedges[j].nlist[2];
04032 jln2=gedges[j].nlist[3];
04033
04034 if (iln1==jfn1 && iln2==jfn2){
04035 gedges[i].next=j;
04036 gedges[j].prev=i;
04037 }
04038 if (ifn1==jln1 && ifn2==jln2){
04039 gedges[i].prev=j;
04040 gedges[j].next=i;
04041 }
04042 if (ifn1==jfn1 && ifn2==jfn2){
04043 gedges[i].prev=j;
04044 gedges[j].prev=i;
04045 }
04046 if (iln1==jln1 && iln2==jln2){
04047 gedges[i].next=j;
04048 gedges[j].next=i;
04049 }
04050
04051 }
04052 }
04053
04054 }
04055
04056
04057
04058
04059
04060
04061
04062
04063 void gtopology::edge_series (FILE *out)
04064 {
04065 long i,j,k,n,c,e,a,m;
04066 long *av;
04067
04068
04069 nser=0;
04070
04071
04072 edgeser = new long [nged];
04073 for (i=0;i<nged;i++){
04074 edgeser[i]=-1;
04075 }
04076
04077 if (nged>0){
04078
04079
04080 c=0;
04081 for (i=0;i<nged;i++){
04082 if (gedges[i].prev==-1 && gedges[i].next==-1){
04083 edgeser[i]=nser;
04084 nser++;
04085 c++;
04086 }
04087 }
04088
04089 a=-2; n=-2;
04090 for (i=0;i<nged;i++){
04091 if (gedges[i].prev==-1 && edgeser[i]==-1){
04092
04093 a=i;
04094 n=gedges[i].next;
04095 edgeser[i]=nser;
04096 break;
04097 }
04098 }
04099 if (a==-2 && n==-2){
04100
04101 for (i=0;i<nged;i++){
04102 if (gedges[i].next==-1 && edgeser[i]==-1){
04103 a=i;
04104 n=gedges[i].prev;
04105 edgeser[i]=nser;
04106
04107 gedges[i].next=gedges[i].prev;
04108 gedges[i].prev=-1;
04109 m=gedges[i].nlist[0];
04110 gedges[i].nlist[0]=gedges[i].nlist[2];
04111 gedges[i].nlist[2]=m;
04112 m=gedges[i].nlist[1];
04113 gedges[i].nlist[1]=gedges[i].nlist[3];
04114 gedges[i].nlist[3]=m;
04115
04116 m=gedges[i].fn;
04117 gedges[i].fn=gedges[i].ln;
04118 gedges[i].ln=m;
04119
04120 break;
04121 }
04122 }
04123 }
04124 if (a==-2 && n==-2){
04125 i=0;
04126 a=i;
04127 n=gedges[i].next;
04128 edgeser[i]=nser;
04129 }
04130 c++;
04131
04132 if (c<nged){
04133 do{
04134 e=gedges[n].next;
04135 if (edgeser[n]==-1){
04136 if (a==e){
04137 gedges[n].next=gedges[n].prev;
04138 gedges[n].prev=a;
04139 m=gedges[n].nlist[0];
04140 gedges[n].nlist[0]=gedges[n].nlist[2];
04141 gedges[n].nlist[2]=m;
04142 m=gedges[n].nlist[1];
04143 gedges[n].nlist[1]=gedges[n].nlist[3];
04144 gedges[n].nlist[3]=m;
04145
04146 m=gedges[n].fn;
04147 gedges[n].fn=gedges[n].ln;
04148 gedges[n].ln=m;
04149
04150 a=n;
04151 n=gedges[n].next;
04152 }
04153 else{
04154 a=n;
04155 n=e;
04156 }
04157 edgeser[a]=nser;
04158 c++;
04159 }else{
04160 n=-1;
04161 }
04162
04163 if (n<0){
04164 j=0;
04165 for (i=0;i<nged;i++){
04166 if (gedges[i].prev==-1 && edgeser[i]==-1){
04167 a=i;
04168 n=gedges[i].next;
04169 nser++;
04170 edgeser[i]=nser;
04171 c++;
04172 j=1;
04173 break;
04174 }
04175 }
04176 if (j==0){
04177 for (i=0;i<nged;i++){
04178 if (gedges[i].next==-1 && edgeser[i]==-1){
04179
04180 a=i;
04181 n=gedges[i].prev;
04182 edgeser[i]=nser;
04183
04184 gedges[i].next=gedges[i].prev;
04185 gedges[i].prev=-1;
04186 m=gedges[i].nlist[0];
04187 gedges[i].nlist[0]=gedges[i].nlist[2];
04188 gedges[i].nlist[2]=m;
04189 m=gedges[i].nlist[1];
04190 gedges[i].nlist[1]=gedges[i].nlist[3];
04191 gedges[i].nlist[3]=m;
04192
04193 m=gedges[i].fn;
04194 gedges[i].fn=gedges[i].ln;
04195 gedges[i].ln=m;
04196
04197 nser++;
04198 edgeser[i]=nser;
04199 c++;
04200 break;
04201 }
04202 }
04203 }
04204 if (j==0){
04205 for (i=0;i<nged;i++){
04206 if (edgeser[i]==-1){
04207 nser++;
04208 edgeser[i]=nser;
04209 c++;
04210 a=i;
04211 n=gedges[i].next;
04212 j=1;
04213 }
04214 if (j==1){
04215 break;
04216 }
04217 }
04218 }
04219 }
04220 }while(c<nged);
04221 nser++;
04222 }
04223 }
04224
04225 fprintf (stdout,"\n\n number of series %ld",nser);
04226 for (i=0;i<nser;i++){
04227 fprintf (out,"%ld\n",edgeser[i]);
04228 }
04229 fprintf (out,"\n");
04230
04231
04232 nedser = new long [nser];
04233 for (i=0;i<nser;i++){
04234 nedser[i]=0;
04235 for (j=0;j<nged;j++){
04236 if (edgeser[j]==i)
04237 nedser[i]++;
04238 }
04239 }
04240
04241
04242
04243
04244
04245 edgelist = new long* [nser];
04246 for (i=0;i<nser;i++){
04247 edgelist[i] = new long [nedser[i]];
04248 }
04249 av = new long [nged];
04250 for (i=0;i<nged;i++){
04251 av[i]=-1;
04252 }
04253
04254 for (i=0;i<nser;i++){
04255
04256
04257 a=-1;
04258 for (j=0;j<nged;j++){
04259 if (av[j]>-1)
04260 continue;
04261 else{
04262 if (gedges[j].prev==-1){
04263 a=j;
04264 break;
04265 }
04266 }
04267 }
04268 if (a==-1){
04269
04270 for (j=0;j<nged;j++){
04271 if (av[j]>-1){
04272 continue;
04273 }
04274 else{
04275 a=j;
04276 break;
04277 }
04278 }
04279 }
04280
04281 k=0;
04282 for (j=0;j<nged;j++){
04283 if (av[a]>-1){
04284
04285 break;
04286 }
04287 edgelist[i][k]=a;
04288 k++;
04289 av[a]=1;
04290 n=gedges[a].next;
04291 a=n;
04292 if (n==-1){
04293
04294 break;
04295 }
04296 }
04297 }
04298
04299
04300 fprintf (out,"\n\n\n kontrola pole edgelist \n");
04301 for (i=0;i<nser;i++){
04302 fprintf (out,"\n serie %5ld ",i);
04303 for (j=0;j<nedser[i];j++){
04304 fprintf (out," %ld %ld %ld\n",edgelist[i][j],gedges[edgelist[i][j]].fn+1,gedges[edgelist[i][j]].ln+1);
04305 }
04306 }
04307 fprintf (out,"\n\n");
04308
04309 delete [] av;
04310 }
04311
04312
04313
04314
04315
04316
04317 void gtopology::edge_elem (FILE *out)
04318 {
04319 long i,j,k,l,m,fn,ln,ef,nne,nc,eid;
04320 ivector nod;
04321
04322
04323 adjacelem (out);
04324
04325
04326 for (i=0;i<nser;i++){
04327
04328 ef=-1;
04329
04330
04331 for (j=0;j<nedser[i];j++){
04332
04333 m=edgelist[i][j];
04334
04335 fn=gedges[m].fn;
04336
04337 ln=gedges[m].ln;
04338
04339 if (ef==-1){
04340
04341
04342 for (k=0;k<ne;k++){
04343
04344 nne = give_original_nne (k);
04345
04346 allocv (nne,nod);
04347
04348 give_original_nodes (k,nod);
04349
04350 nc=0;
04351
04352 for (l=0;l<nne;l++){
04353 if (nod[l]==fn || nod[l]==ln)
04354 nc++;
04355 }
04356 destrv (nod);
04357
04358 if (nc==2){
04359 gedges[m].re=k;
04360 ef=k;
04361 break;
04362 }
04363 }
04364 }
04365 else{
04366
04367
04368
04369
04370 for (k=0;k<nadjelel[ef];k++){
04371
04372 eid=adjelel[ef][k];
04373
04374 nne = give_original_nne (eid);
04375
04376 allocv (nne,nod);
04377
04378 give_original_nodes (eid,nod);
04379
04380 nc=0;
04381
04382 for (l=0;l<nne;l++){
04383 if (nod[l]==fn || nod[l]==ln)
04384 nc++;
04385 }
04386 destrv (nod);
04387
04388 if (nc==2){
04389 gedges[m].re=eid;
04390 ef=eid;
04391 break;
04392 }
04393
04394 }
04395 }
04396
04397 }
04398 }
04399 }
04400
04401
04402
04403
04404
04405
04406
04407 void gtopology::edge_allelem (FILE *out)
04408 {
04409 long i,j,k,l,m,fn,ln,ef,nne,nc,eid;
04410 ivector nod;
04411
04412 for (i=0;i<nged;i++){
04413 gedges[i].adjel = new long [2];
04414 gedges[i].adjel[0]=gedges[i].re;
04415 }
04416
04417
04418 for (i=0;i<nser;i++){
04419
04420 ef=-1;
04421
04422
04423 for (j=0;j<nedser[i];j++){
04424
04425 m=edgelist[i][j];
04426
04427 fn=gedges[m].nlist[1];
04428
04429 ln=gedges[m].nlist[3];
04430
04431 if (ef==-1){
04432
04433
04434 for (k=0;k<ne;k++){
04435
04436 nne = give_original_nne (k);
04437
04438 allocv (nne,nod);
04439
04440 give_original_nodes (k,nod);
04441
04442 nc=0;
04443
04444 for (l=0;l<nne;l++){
04445 if (nod[l]==fn || nod[l]==ln)
04446 nc++;
04447 }
04448 destrv (nod);
04449
04450 if (nc==2 && gedges[m].adjel[0]!=k){
04451 gedges[m].adjel[1]=k;
04452 ef=k;
04453 break;
04454 }
04455 }
04456 }
04457 else{
04458
04459
04460
04461
04462 for (k=0;k<nadjelel[ef];k++){
04463
04464 eid=adjelel[ef][k];
04465
04466 nne = give_original_nne (eid);
04467
04468 allocv (nne,nod);
04469
04470 give_original_nodes (eid,nod);
04471
04472 nc=0;
04473
04474 for (l=0;l<nne;l++){
04475 if (nod[l]==fn || nod[l]==ln)
04476 nc++;
04477 }
04478 destrv (nod);
04479
04480 if (nc==2 && gedges[m].adjel[0]!=eid){
04481 gedges[m].adjel[1]=eid;
04482 ef=eid;
04483 break;
04484 }
04485
04486 }
04487 }
04488
04489 }
04490 }
04491
04492 }
04493
04494
04495
04496
04497
04498
04499
04500
04501 void gtopology::normvectorient ()
04502 {
04503 long i,eid,nne;
04504 ivector nod;
04505 vector x,y;
04506
04507 for (i=0;i<nged;i++){
04508
04509 eid=gedges[i].re;
04510
04511 nne = give_original_nne (eid);
04512
04513 allocv (nne,x);
04514 allocv (nne,y);
04515 allocv (nne,nod);
04516
04517
04518 give_node_coord2d (x,y,eid);
04519
04520 give_original_nodes (eid,nod);
04521
04522
04523 gedges[i].check_normal (x,y,nod);
04524
04525
04526 destrv (nod);
04527 destrv (x);
04528 destrv (y);
04529 }
04530
04531 }
04532
04533
04534
04535
04536
04537
04538
04539
04540
04541
04542
04543
04544
04545
04546
04547
04548
04549
04550
04551
04552
04553
04554
04555
04556
04557
04558
04559
04560
04561
04562
04563
04564
04565
04566
04567
04568
04569
04570
04571
04572
04573
04574
04575
04576
04577
04578
04579
04580
04581
04582
04583
04584
04585
04586
04587
04588
04589
04590
04591
04592
04593
04594
04595
04596
04597
04598
04599
04600
04601
04602
04603
04604
04605
04606
04607
04608
04609
04610
04611
04612
04613
04614
04615
04616
04617
04618
04619
04620
04621
04622
04623
04624
04625
04626
04627 void gtopology::edges_auxprint (FILE *out)
04628 {
04629 long i,j,a;
04630
04631 fprintf (out,"\n\n\n kontrola pole nlist na jednotlivych hranach \n");
04632 for (i=0;i<nged;i++){
04633 fprintf (out,"\n %6ld %6ld %6ld %6ld",gedges[i].nlist[0]+1,gedges[i].nlist[1]+1,gedges[i].nlist[2]+1,gedges[i].nlist[3]+1);
04634 }
04635
04636
04637
04638 fprintf (out,"\n\n\n kontrola usporadni hran \n");
04639 for (i=0;i<nged;i++){
04640 fprintf (out,"\n %6ld %6ld",gedges[i].nlist[0]+1,gedges[i].nlist[1]+1);
04641 }
04642 i=nged-1;
04643 if (i>=0)
04644 fprintf (out,"\n %6ld %6ld",gedges[i].nlist[2]+1,gedges[i].nlist[3]+1);
04645
04646
04647
04648 fprintf (out,"\n\n\n kontrola prev a next \n");
04649 for (i=0;i<nged;i++){
04650 fprintf (out,"\n hrana %6ld prev %6ld next %6ld",i,gedges[i].prev+1,gedges[i].next+1);
04651 }
04652
04653
04654 fprintf (out,"\n\n\n kontrola souslednosti \n");
04655 for (i=0;i<nser;i++){
04656 fprintf (out,"\n\n serie %5ld",i);
04657 for (j=0;j<nedser[i];j++){
04658 a=edgelist[i][j];
04659 fprintf (out,"\n %5ld hrana %5ld %5ld %5ld %lf %lf %lf %lf %lf %lf",
04660 j,a,gedges[a].fn,gedges[a].ln,
04661 gedges[a].dv[0],gedges[a].dv[1],gedges[a].dv[2],
04662 gedges[a].nv[0],gedges[a].nv[1],gedges[a].nv[2]);
04663 }
04664 }
04665 fprintf (out,"\n\n");
04666
04667 fprintf (out,"\n\n\n kontrola prilehlych prvku \n");
04668 for (i=0;i<nser;i++){
04669 fprintf (out,"\n\n serie %5ld",i);
04670 for (j=0;j<nedser[i];j++){
04671 a=edgelist[i][j];
04672 fprintf (out,"\n hrana %5ld %6ld %6ld",j,gedges[a].adjel[0],gedges[a].adjel[1]);
04673 }
04674 }
04675 fprintf (out,"\n\n");
04676
04677
04678
04679
04680
04681
04682
04683
04684
04685
04686
04687
04688
04689
04690
04691
04692
04693
04694
04695
04696
04697
04698
04699
04700
04701
04702
04703
04704
04705
04706
04707
04708
04709
04710
04711
04712
04713
04714
04715
04716
04717
04718
04719
04720
04721
04722
04723
04724
04725
04726
04727
04728
04729
04730
04731
04732
04733
04734
04735
04736
04737
04738
04739
04740
04741
04742
04743
04744
04745
04746
04747
04748
04749
04750
04751
04752
04753
04754
04755
04756
04757
04758
04759
04760
04761
04762
04763
04764
04765
04766
04767
04768
04769
04770
04771
04772
04773
04774
04775
04776
04777
04778
04779
04780
04781
04782
04783
04784
04785
04786
04787
04788
04789
04790
04791
04792
04793
04794
04795
04796
04797
04798
04799
04800
04801
04802
04803
04804
04805
04806
04807
04808
04809
04810
04811
04812
04813
04814
04815
04816
04817
04818
04819
04820
04821
04822
04823
04824
04825
04826
04827
04828
04829
04830
04831
04832
04833
04834
04835
04836
04837
04838
04839
04840
04841
04842
04843
04844
04845
04846
04847
04848
04849
04850
04851
04852
04853
04854
04855
04856
04857
04858
04859
04860
04861
04862
04863
04864
04865
04866
04867
04868
04869
04870
04871
04872
04873
04874
04875
04876
04877
04878
04879
04880
04881
04882
04883
04884
04885
04886
04887
04888
04889
04890
04891
04892
04893 }
04894
04895
04896
04897
04898
04899
04900
04901
04902
04903
04904
04905
04906
04907
04908
04909
04910
04911
04912
04913
04914
04915
04916
04917
04918
04919
04920
04921
04922
04923
04924
04925
04926
04927
04928
04929
04930
04931
04932
04933
04934
04935
04936
04937
04938
04939
04940
04941
04942
04943
04944
04945
04946
04947
04948
04949
04950
04951
04952
04953
04954
04955
04956 void gtopology::alloc_growstr()
04957 {
04958 long i,j,ndofn;
04959
04960 for(i=0; i<nn; i++)
04961 {
04962 if (gnodes[i].tgf == NULL)
04963 {
04964
04965 ndofn = gnodes[i].give_ndofn();
04966 if (ndofn > 0)
04967 {
04968 gnodes[i].tgf = new long[ndofn];
04969
04970
04971 for(j=0; j<ndofn; j++)
04972 gnodes[i].tgf[j] = -1;
04973 }
04974 }
04975 }
04976 }
04977
04978
04979
04980
04981
04982
04983
04984
04985
04986 void gtopology::alloc_edge_cn ()
04987 {
04988 long i,nid,ndofnfn,ndofnln;
04989
04990 for (i=0;i<nged;i++){
04991
04992 nid=gedges[i].fn;
04993
04994 ndofnfn=give_ndofn (nid);
04995
04996 nid=gedges[i].ln;
04997
04998 ndofnln=give_ndofn (nid);
04999
05000
05001 gedges[i].alloc_cn (ndofnfn,ndofnln);
05002 }
05003
05004 }
05005
05006
05007
05008
05009
05010
05011
05012
05013
05014 void gtopology::edge_localization (gmatrix *gm)
05015 {
05016 long i,nm,*ncn1,*ncn2,*mcn;
05017
05018 for (i=0;i<nged;i++){
05019
05020 nm=gedges[i].ndofnf;
05021 ncn1 = new long [nm];
05022 ncn2 = new long [nm];
05023 mcn = new long [nm];
05024
05025
05026 give_edge_code_numbers (i,1,ncn1,ncn2,mcn);
05027
05028
05029 gm->mult_localize (nm,ncn1,ncn2,mcn);
05030
05031 delete [] mcn;
05032 delete [] ncn2;
05033 delete [] ncn1;
05034
05035
05036
05037 nm=gedges[i].ndofnl;
05038 ncn1 = new long [nm];
05039 ncn2 = new long [nm];
05040 mcn = new long [nm];
05041
05042
05043 give_edge_code_numbers (i,2,ncn1,ncn2,mcn);
05044
05045
05046 gm->mult_localize (nm,ncn1,ncn2,mcn);
05047
05048 delete [] mcn;
05049 delete [] ncn2;
05050 delete [] ncn1;
05051 }
05052
05053 }
05054
05055
05056
05057
05058
05059
05060
05061
05062
05063
05064
05065
05066
05067
05068
05069
05070
05071
05072
05073
05074
05075
05076
05077
05078
05079
05080
05081
05082
05083
05084
05085
05086 void gtopology::auto_subdom_detection (FILE *out)
05087 {
05088
05089
05090
05091
05092
05093 end_nodes (out);
05094
05095
05096 enodes_allelem (out);
05097
05098 endnodes_auxprint (out);
05099
05100
05101
05102
05103
05104
05105
05106 if (edtype==jumps){
05107
05108 adjacnodes (out);
05109
05110
05111 edges (out);
05112 }
05113 if (edtype==mater){
05114
05115 adjacelem (out);
05116
05117
05118 long matrad=251;
05119
05120 edges_mat (matrad,out);
05121 }
05122
05123
05124 edgenode_sorting ();
05125
05126 prev_next_edges ();
05127
05128 edge_series (out);
05129
05130 edge_elem (out);
05131
05132
05133
05134 edge_allelem (out);
05135
05136
05137 edge_dirvect ();
05138
05139 edge_normvect ();
05140
05141 normvectorient ();
05142
05143
05144 if (edtype==jumps){
05145
05146 domdecomp ();
05147
05148
05149 create_ltg (out);
05150 }
05151
05152 edges_auxprint (out);
05153
05154 if (edtype==mater){
05155 cngen=1;
05156 }
05157 if (edtype==jumps){
05158
05159 cngen=3;
05160 }
05161
05162 }
05163
05164
05165
05166
05167
05168
05169
05170 long gtopology::codenum_multip ()
05171 {
05172 long i,j,k,l,ii,nid,nm,ndofn;
05173 long *cn;
05174
05175
05176 nm=ndof;
05177
05178
05179
05180
05181
05182
05183 alloc_endnode_cn ();
05184
05185 for (i=0;i<nen;i++){
05186
05187 nid=endnodes[i].fn;
05188
05189 ndofn=give_ndofn (nid);
05190 cn = new long [ndofn];
05191 for (l=0;l<ndofn;l++){
05192 ii=give_dof (nid,l);
05193 if (ii!=0){
05194 cn[l]=nm;
05195 nm++;
05196 }
05197 else
05198 cn[l]=0;
05199 }
05200 for (l=0;l<ndofn;l++){
05201 endnodes[i].cnm[l]=cn[l];
05202 }
05203 delete [] cn;
05204 }
05205
05206
05207
05208
05209
05210
05211 alloc_edge_cn ();
05212
05213 for (i=0;i<nser;i++){
05214
05215 for (j=0;j<nedser[i];j++){
05216
05217
05218
05219 k=edgelist[i][j];
05220
05221
05222
05223
05224 if (j==0){
05225
05226 nid=gedges[k].fn;
05227
05228 ndofn=give_ndofn (nid);
05229 cn = new long [ndofn];
05230 for (l=0;l<ndofn;l++){
05231 ii=give_dof (nid,l);
05232 if (ii!=0){
05233 cn[l]=nm;
05234 nm++;
05235 }
05236 else
05237 cn[l]=0;
05238 }
05239 }
05240
05241 for (l=0;l<ndofn;l++){
05242 gedges[k].cnfn[l]=cn[l];
05243 }
05244 delete [] cn;
05245
05246
05247
05248
05249
05250 nid=gedges[k].ln;
05251
05252 ndofn=give_ndofn (nid);
05253 cn = new long [ndofn];
05254 for (l=0;l<ndofn;l++){
05255 ii=give_dof (nid,l);
05256 if (ii!=0){
05257 cn[l]=nm;
05258 nm++;
05259 }
05260 else
05261 cn[l]=0;
05262 }
05263
05264 for (l=0;l<ndofn;l++){
05265 gedges[k].cnln[l]=cn[l];
05266 }
05267
05268 if (j==nedser[i]-1)
05269 delete [] cn;
05270 }
05271 }
05272
05273 return nm;
05274 }
05275
05276
05277
05278
05279
05280
05281
05282
05283
05284
05285
05286 void gtopology::mult_localization (gmatrix *gm)
05287 {
05288
05289 endnodes_localization (gm);
05290
05291 edge_localization (gm);
05292 }
05293
05294
05295
05296
05297
05298
05299
05300 void gtopology::domdecomp ()
05301 {
05302 long i,j,ii,jj,nvn,nan,nnf,nnnf,anode;
05303 long *aux,*oldlist,*newlist;
05304
05305 aux = new long [nn];
05306 for (i=0;i<nn;i++){
05307 aux[i]=-1;
05308 }
05309
05310
05311 ns=0;
05312
05313 nvn=0;
05314 do{
05315 for (i=0;i<nn;i++){
05316 if (aux[i]==-1)
05317 anode=i;
05318 }
05319
05320 oldlist=new long[1];
05321 oldlist[0]=anode;
05322 nnf=1;
05323
05324 do{
05325
05326 nnnf=0;
05327
05328
05329 for (i=0;i<nnf;i++){
05330
05331 ii=oldlist[i];
05332
05333 nan=nadjnodnod[ii];
05334
05335 for (j=0;j<nan;j++){
05336
05337 jj=adjnodnod[ii][j];
05338 if (aux[jj]==-1){
05339
05340 nnnf++;
05341 }
05342 }
05343 }
05344
05345 newlist=new long [nnnf];
05346
05347 nnnf=0;
05348
05349
05350 for (i=0;i<nnf;i++){
05351
05352 ii=oldlist[i];
05353
05354 nan=nadjnodnod[ii];
05355
05356 for (j=0;j<nan;j++){
05357
05358 jj=adjnodnod[ii][j];
05359 if (aux[jj]==-1){
05360
05361 newlist[nnnf]=jj;
05362 nnnf++;
05363 aux[jj]=ns;
05364 nvn++;
05365 }
05366 }
05367 }
05368 delete [] oldlist;
05369 oldlist=new long [nnnf];
05370 for (i=0;i<nnnf;i++){
05371 oldlist[i]=newlist[i];
05372 }
05373 delete [] newlist;
05374 nnf=nnnf;
05375 }while(nnf>0);
05376
05377 delete [] oldlist;
05378
05379 ns++;
05380 }while(nvn<nn);
05381
05382
05383
05384 if (stop==NULL)
05385 stop = new seqtop (ns,all_nodes);
05386
05387
05388 if (stop->nnsd==NULL)
05389 stop->nnsd = new long [ns];
05390 for (i=0;i<ns;i++){
05391 stop->nnsd[i]=0;
05392 }
05393
05394 for (i=0;i<nn;i++){
05395 stop->nnsd[aux[i]]++;
05396 }
05397
05398 delete [] aux;
05399 }
05400
05401
05402
05403
05404
05405
05406
05407
05408
05409 void gtopology::create_ltg (FILE *out)
05410 {
05411 long i,j,k,nin;
05412 long *aux,*v;
05413
05414 aux = new long [nn];
05415 for (i=0;i<nn;i++){
05416 aux[i]=nn;
05417 }
05418
05419
05420
05421
05422
05423
05424
05425
05426 for (i=0;i<nen;i++){
05427 k=endnodes[i].fn;
05428 j=endnodes[i].ln;
05429 if (j<k){
05430 if (aux[j]>j)
05431 aux[j]=j;
05432 if (aux[k]>j)
05433 aux[k]=j;
05434
05435 if (aux[j]<aux[k])
05436 aux[k]=aux[j];
05437 else
05438 aux[j]=aux[k];
05439 }
05440 else{
05441 if (aux[j]>k)
05442 aux[j]=k;
05443 if (aux[k]>k)
05444 aux[k]=k;
05445
05446 if (aux[j]<aux[k])
05447 aux[k]=aux[j];
05448 else
05449 aux[j]=aux[k];
05450 }
05451
05452 }
05453
05454
05455 for (i=0;i<nged;i++){
05456 j=gedges[i].nlist[0];
05457 k=gedges[i].nlist[1];
05458 if (j<k){
05459 if (aux[j]>j)
05460 aux[j]=j;
05461 if (aux[k]>j)
05462 aux[k]=j;
05463
05464 if (aux[j]<aux[k])
05465 aux[k]=aux[j];
05466 else
05467 aux[j]=aux[k];
05468 }
05469 else{
05470 if (aux[j]>k)
05471 aux[j]=k;
05472 if (aux[k]>k)
05473 aux[k]=k;
05474
05475 if (aux[j]<aux[k])
05476 aux[k]=aux[j];
05477 else
05478 aux[j]=aux[k];
05479 }
05480
05481 j=gedges[i].nlist[2];
05482 k=gedges[i].nlist[3];
05483 if (j<k){
05484 if (aux[j]>j)
05485 aux[j]=j;
05486 if (aux[k]>j)
05487 aux[k]=j;
05488
05489 if (aux[j]<aux[k])
05490 aux[k]=aux[j];
05491 else
05492 aux[j]=aux[k];
05493 }
05494 else{
05495 if (aux[j]>k)
05496 aux[j]=k;
05497 if (aux[k]>k)
05498 aux[k]=k;
05499
05500 if (aux[j]<aux[k])
05501 aux[k]=aux[j];
05502 else
05503 aux[j]=aux[k];
05504 }
05505
05506 }
05507
05508
05509
05510 v = new long [nn];
05511 for (i=0;i<nn;i++){
05512 v[i]=-1;
05513 }
05514
05515
05516 for (i=0;i<nn;i++){
05517 j=aux[i];
05518 if (j!=nn)
05519 v[j]++;
05520 }
05521
05522
05523
05524 nin=0;
05525 for (i=0;i<nn;i++){
05526 if (v[i]>-1){
05527 v[i]=nin;
05528 nin++;
05529 }
05530 }
05531
05532 for (i=0;i<nn;i++){
05533 j=aux[i];
05534 if (j!=nn)
05535 aux[i]=v[j];
05536 else
05537 aux[i]=-1;
05538 }
05539
05540
05541
05542
05543
05544 if (stop->ltg==NULL){
05545 stop->ltg = new long* [stop->ns];
05546 for (i=0;i<stop->ns;i++){
05547 stop->ltg[i] = new long [stop->nnsd[i]];
05548 }
05549 }
05550
05551 k=0;
05552 for (i=0;i<stop->ns;i++){
05553 for (j=0;j<stop->nnsd[i];j++){
05554 stop->ltg[i][j]=aux[k];
05555 k++;
05556 }
05557 }
05558
05559
05560 fprintf (out,"\n\n\n funkce seqtop::create_ltg pole ltg");
05561 for (i=0;i<stop->ns;i++){
05562 fprintf (out,"\n");
05563 for (j=0;j<stop->nnsd[i];j++){
05564 fprintf (out,"%ld ",stop->ltg[i][j]);
05565 }
05566 }
05567 fprintf (out,"\n");
05568
05569
05570 delete [] aux;
05571 delete [] v;
05572
05573 }
05574
05575
05576
05577
05578
05579
05580
05581
05582
05583
05584
05585
05586 void gtopology::read_seq_top (XFILE *in)
05587 {
05588 meshdescription md;
05589
05590
05591
05592
05593
05594
05595
05596
05597 xfscanf (in,"%k%m","meshdescript",&meshdescription_kwdset,(int*)&md);
05598
05599
05600 xfscanf (in,"%ld",&ns);
05601
05602
05603 stop = new seqtop (ns,md);
05604
05605 if (md==all_nodes || md==bound_nodes || md==neg_bound_nodes || md==metis){
05606
05607
05608
05609
05610
05611 stop->read_nnsd (in);
05612
05613
05614 stop->read_ltg (in);
05615 }
05616 if (md==metis_elem){
05617
05618
05619
05620
05621
05622
05623
05624
05625
05626
05627 stop->read_eldom (ne,in);
05628 }
05629
05630 }
05631
05632
05633
05634
05635
05636
05637
05638
05639
05640
05641
05642
05643
05644
05645
05646
05647
05648
05649
05650
05651
05652
05653
05654
05655
05656
05657
05658
05659
05660
05661
05662
05663
05664
05665
05666
05667
05668
05669
05670
05671
05672
05673
05674
05675
05676
05677
05678
05679
05680
05681
05682
05683
05684
05685
05686
05687
05688
05689
05690
05691
05692
05693
05694
05695
05696
05697
05698
05699
05700
05701
05702
05703
05704
05705
05706
05707
05708
05709
05710
05711
05712
05713
05714
05715
05716
05717
05718
05719
05720
05721
05722
05723
05724
05725
05726
05727
05728
05729
05730
05731
05732
05733
05734
05735
05736
05737
05738
05739
05740
05741
05742
05743
05744
05745
05746
05747
05748
05749
05750
05751
05752
05753
05754
05755
05756
05757
05758
05759
05760
05761
05762
05763
05764
05765
05766
05767
05768
05769
05770
05771
05772
05773
05774
05775
05776
05777
05778
05779
05780
05781
05782
05783
05784
05785
05786
05787
05788
05789
05790
05791
05792
05793
05794
05795
05796
05797
05798
05799
05800
05801
05802
05803
05804
05805
05806
05807
05808
05809
05810
05811
05812
05813
05814
05815 void gtopology::cuthill_mckee_renumb (FILE *out)
05816 {
05817 long i,j,k;
05818 long minnode,n,a,b,m;
05819 int startnode;
05820 long min;
05821 long *cmk;
05822 bool stop;
05823
05824
05825
05826 adjacnodes_edge(out);
05827
05828
05829 cmk = new long [nn];
05830 for(i = 0; i<nn; i++) cmk[i] = -1;
05831
05832
05833
05834 min = LONG_MAX;
05835 for(i = 0; i<nn; i++){
05836 if(min >= nadjacnodesedge[i]){
05837 min = nadjacnodesedge[i];
05838 minnode = i;
05839 }
05840 }
05841 fprintf (stdout,"\n\n\nminnode = %ld\n\n",minnode);
05842
05843
05844
05845 cmk[minnode] = 0;
05846
05847
05848 for(i = 0; i < nn; i++){
05849
05850 for (m = nadjacnodesedge[i] - 1; m >= 0; m--){
05851 for (j = 1; j <= m; j++){
05852 a = adjacnodesedge[i][j];
05853 b = adjacnodesedge[i][j-1];
05854 if (nadjacnodesedge[b] > nadjacnodesedge[a]){
05855 n = adjacnodesedge[i][j-1];
05856 adjacnodesedge[i][j-1] = adjacnodesedge[i][j];
05857 adjacnodesedge[i][j] = n;
05858 }
05859 }
05860 }
05861 }
05862 for(i = 0; i < nn; i++){
05863 fprintf(stdout,"node %ld # neighbours %ld:",i,nadjacnodesedge[i]);
05864 for(j= 0; j < nadjacnodesedge[i]; j++){
05865 fprintf(stdout," %ld",adjacnodesedge[i][j]);
05866 }
05867 fprintf(stdout,"\n");
05868 }
05869
05870 startnode = 0;
05871 stop = 0;
05872 m = 1;
05873 stop = false;
05874 for(k = 0; k < nn; k++){
05875 for(i = 0; i < nn; i++){
05876 if(cmk[i] == startnode){
05877
05878 for(j = 0; j < nadjacnodesedge[i];j++){
05879 if(cmk[adjacnodesedge[i][j]] == -1){
05880
05881 cmk[adjacnodesedge[i][j]] = m;
05882
05883 m++;
05884 }
05885 if(m == nn){
05886 stop = true;
05887 break;
05888 }
05889 }
05890 break;
05891 }
05892 }
05893 if(stop == true){
05894 break;
05895 }
05896 startnode++;
05897 }
05898
05899
05900
05901
05902
05903 for(i = 0; i < nn; i++) delete[] adjacnodesedge[i];
05904 delete []adjacnodesedge;
05905 delete []nadjacnodesedge;
05906
05907
05908
05909
05910 if (nodren == cuthill_mckee){
05911 for (i=0;i<nn;i++){
05912 ordering[cmk[i]]=i;
05913 }
05914 }
05915
05916 if (nodren == rev_cuthill_mckee){
05917 for (i=0;i<nn;i++){
05918 ordering[cmk[i]]=i;
05919 }
05920 j=nn-1;
05921 for (i=0;i<nn;i++){
05922 cmk[j]=ordering[i];
05923 j--;
05924 }
05925 for (i=0;i<nn;i++){
05926 ordering[i]=cmk[i];
05927 }
05928
05929 }
05930
05931
05932
05933
05934
05935
05936
05937 }
05938
05939
05940
05941
05942
05943 void gtopology::writePriorityQueue(int *fronta){
05944 for (int i = 0; i < (nn); i++){
05945
05946 }
05947 }
05948
05949 void gtopology::shell_sort( int *array, int arrayLength)
05950 {
05951 int flag = 1, d = arrayLength, i, temp;
05952 while( flag || (d>1))
05953 {
05954 flag = 0;
05955 d = (d+1) / 2;
05956 for (i = 0; i < (arrayLength - d); i++)
05957 {
05958 if (nadjnodnod[array[i + d]] < nadjnodnod[array[i]])
05959 {
05960 temp = array[i + d];
05961 array[i + d] = array[i];
05962 array[i] = temp;
05963 flag = 1;
05964 }
05965 }
05966 }
05967 return;
05968 }
05969
05970
05971
05972
05973
05974
05975
05976
05977
05978
05979
05980
05981
05982 void gtopology::shell_sort_x(ivector &x_ord)
05983 {
05984 long flag = 1, d = nn, i, temp;
05985
05986 reallocv(nn, x_ord);
05987 for(i=0; i<nn; i++)
05988 x_ord[i] = i;
05989
05990 while(flag || (d>1))
05991 {
05992 flag = 0;
05993 d = (d+1)/2;
05994 for (i=0; i<(nn-d); i++)
05995 {
05996 if (gnodes[i+d].x < gnodes[i].x)
05997 {
05998 temp = x_ord[i+d];
05999 x_ord[i+d] = x_ord[i];
06000 x_ord[i] = temp;
06001 flag = 1;
06002 }
06003 }
06004 }
06005 return;
06006 }
06007
06008
06009
06010
06011
06012
06013
06014
06015
06016
06017
06018
06019
06020 void gtopology::shell_sort_y(ivector &y_ord)
06021 {
06022 long flag = 1, d = nn, i, temp;
06023
06024 reallocv(nn, y_ord);
06025 for(i=0; i<nn; i++)
06026 y_ord[i] = i;
06027
06028 while(flag || (d>1))
06029 {
06030 flag = 0;
06031 d = (d+1)/2;
06032 for (i=0; i<(nn-d); i++)
06033 {
06034 if (gnodes[i+d].y < gnodes[i].y)
06035 {
06036 temp = y_ord[i+d];
06037 y_ord[i+d] = y_ord[i];
06038 y_ord[i] = temp;
06039 flag = 1;
06040 }
06041 }
06042 }
06043 return;
06044 }
06045
06046
06047
06048
06049
06050
06051
06052
06053
06054
06055
06056
06057
06058 void gtopology::shell_sort_z(ivector &z_ord)
06059 {
06060 long flag = 1, d = nn, i, temp;
06061
06062 reallocv(nn, z_ord);
06063 for(i=0; i<nn; i++)
06064 z_ord[i] = i;
06065
06066 while(flag || (d>1))
06067 {
06068 flag = 0;
06069 d = (d+1)/2;
06070 for (i=0; i<(nn-d); i++)
06071 {
06072 if (gnodes[i+d].z < gnodes[i].z)
06073 {
06074 temp = z_ord[i+d];
06075 z_ord[i+d] = z_ord[i];
06076 z_ord[i] = temp;
06077 flag = 1;
06078 }
06079 }
06080 }
06081 return;
06082 }
06083
06084
06085
06086
06087 void gtopology::lastLevelOfRootedStructure(int start, int *velikost, int *maxSirka, int *hloubka, int *vzdalenost, int **lastLevel){
06088 int *aktUzlyStart;
06089 int i,j,k;
06090 aktUzlyStart = new int[nn];
06091 for(i = 0; i< nn; i++){
06092 aktUzlyStart[i] = 0;
06093 }
06094 int pocetObjevenych = 0;
06095 int *hranice;
06096 int velikostHranice;
06097 int aktualniUzel,pomUzel,posl;
06098
06099 *maxSirka = -1;
06100 *hloubka = 0;
06101 aktUzlyStart[start]=1;
06102 velikostHranice = 1;
06103 hranice = new int[1];
06104 hranice[0] = start;
06105 int *pomocnePole;
06106 pomocnePole = new int[nn];
06107 do{
06108
06109
06110
06111 for( j=0; j<nn; j++){
06112 pomocnePole[j] = 0;
06113
06114 }
06115 posl =0;
06116 pocetObjevenych = 0;
06117
06118 for( j=0; j<velikostHranice; j++){
06119
06120 aktualniUzel = hranice[j];
06121
06122
06123 for(i = 0; i< nadjnodnod[aktualniUzel]; i++){
06124 pomUzel = adjnodnod[aktualniUzel][i];
06125 vzdalenost[aktualniUzel] = *hloubka;
06126
06127 if(aktUzlyStart[pomUzel] == 0){
06128
06129
06130 aktUzlyStart[pomUzel] = 1;
06131 pomocnePole[posl] = pomUzel;
06132 pocetObjevenych++;
06133 posl++;
06134
06135 }
06136
06137 }
06138 }
06139 if(pocetObjevenych > *maxSirka){
06140 *maxSirka = pocetObjevenych;
06141 }
06142 (*hloubka)++;
06143
06144
06145 if(pocetObjevenych != 0){
06146 delete[] hranice;
06147 hranice = new int[posl + 1];
06148 for( k=0; k<posl; k++){
06149 hranice[k]= pomocnePole[k];
06150 }
06151 velikostHranice = posl ;
06152
06153 }
06154 }while(pocetObjevenych != 0);
06155
06156
06157
06158
06159
06160 if(*lastLevel != NULL){
06161 delete [] *lastLevel;
06162 }
06163 *lastLevel = hranice;
06164 *velikost = velikostHranice;
06165
06166 }
06167
06168 void gtopology::generateInitialList(int start, int **seznam, int * velikostSeznamu, int *hloubkaSeznamu, int *maxSirka){
06169
06170
06171 int sirka, hloubka,velikost;
06172 int *lastLevel=NULL;
06173 int *vzdalenosti;
06174
06175
06176 vzdalenosti = new int[nn];
06177
06178
06179
06180 lastLevelOfRootedStructure(start,&velikost,&sirka , &hloubka, vzdalenosti , &lastLevel);
06181
06182
06183
06184
06185 shell_sort(lastLevel,velikost);
06186 removeRepetitiousDegrees(&lastLevel, &velikost);
06187
06188 *velikostSeznamu = velikost;
06189 *hloubkaSeznamu = hloubka;
06190 *maxSirka = sirka;
06191 *seznam = lastLevel ;
06192 }
06193
06194
06195
06196 int gtopology::findMinimumDegree(){
06197 int minUzel,minimum = 99999999;
06198 for(int i = 0; i< nn; i++){
06199
06200
06201
06202 if(nadjnodnod[i]< minimum){
06203 minUzel = i;
06204 minimum = nadjnodnod[i];
06205 }
06206 }
06207 return minUzel;
06208 }
06209
06210
06211
06212 void gtopology::removeRepetitiousDegrees(int **pole, int *velikost){
06213 int aktStupen= -1,k,posl;
06214 int *pomocnePole;
06215 int *pracovniPole;
06216 pracovniPole= *pole;
06217 pomocnePole = new int[nn];
06218
06219 for( k=0; k<nn; k++){
06220 pomocnePole[k] = 0;
06221
06222 }
06223 posl =0;
06224
06225
06226 for( k=0; k<*velikost; k++){
06227 if(aktStupen != nadjnodnod[pracovniPole[k]]){
06228
06229 aktStupen = nadjnodnod[pracovniPole[k]];
06230
06231 pomocnePole[posl] = pracovniPole[k];
06232
06233 posl++;
06234 }
06235 }
06236 delete[] pracovniPole;
06237 pracovniPole = new int[posl + 1];
06238 for( k=0; k<posl; k++){
06239 pracovniPole[k]= pomocnePole[k];
06240 }
06241 *velikost = posl ;
06242 *pole = pracovniPole;
06243
06244 }
06245
06246 void gtopology::findPseudoPheripheral(int *startUzel, int *cilUzel){
06247 int *pracovniSeznam;
06248 int velikostSeznamu;
06249 int hloubkaSeznamu;
06250 int maxSirkaSeznamu;
06251 int start,cil ;
06252 int aktualniSirka = 99999999;
06253 int zpracSirka,zpracHloubka,zpracVelikost,*zpracVzdalenosti,*zpracLastLevel;
06254 bool nalezenCil = false;
06255 bool novyStart = false;
06256 zpracLastLevel = NULL;
06257 start = findMinimumDegree();
06258
06259
06260 zpracVzdalenosti = new int[nn];
06261 while(!((nalezenCil)&(!novyStart))){
06262
06263
06264
06265 generateInitialList(start, &pracovniSeznam, &velikostSeznamu, &hloubkaSeznamu, &maxSirkaSeznamu);
06266
06267
06268
06269
06270 novyStart = false;
06271 nalezenCil = false;
06272
06273 for(int k=0; k<velikostSeznamu; k++){
06274
06275
06276
06277 lastLevelOfRootedStructure(pracovniSeznam[k], &zpracVelikost, &zpracSirka, &zpracHloubka, zpracVzdalenosti, &zpracLastLevel) ;
06278
06279
06280
06281
06282 if((zpracHloubka > hloubkaSeznamu) && (zpracSirka < aktualniSirka)){
06283
06284
06285 start = pracovniSeznam[k];
06286 novyStart = true;
06287
06288 delete[] pracovniSeznam;
06289 break;
06290
06291 }
06292 else if(zpracSirka < aktualniSirka){
06293
06294 cil = pracovniSeznam[k];
06295
06296 aktualniSirka = zpracSirka;
06297 nalezenCil = true;
06298 break;
06299 }
06300
06301
06302 }
06303
06304 }
06305 *startUzel = start;
06306 *cilUzel = cil;
06307
06308 }
06309
06310 void gtopology::sloan_renumb (FILE *out)
06311 {
06312 adjacnodes (out);
06313
06314 int zpracSirka,zpracHloubka,zpracVelikost,*zpracLastLevel;
06315
06316 zpracLastLevel = NULL;
06317 int start,cil;
06318 int oznacenych=-1;
06319 int *vzdalenosti, *priority, *statusy, *priorFronta,*noveCislovani;
06320 int vrcholFronty;
06321 int zpracUzel;
06322 int maxPriorita;
06323 int pomUzel,pomActive, poziceUzlu;
06324 int j,k;
06325
06326
06327 findPseudoPheripheral(&start, &cil);
06328
06329
06330
06331
06332
06333
06334
06335 vzdalenosti = new int[nn];
06336 priority = new int[nn];
06337 statusy = new int[nn];
06338 priorFronta = new int[nn];
06339 noveCislovani = new int[nn];
06340 int W1 = 1, W2 = 2;
06341
06342 lastLevelOfRootedStructure(cil, &zpracVelikost, &zpracSirka, &zpracHloubka, vzdalenosti, &zpracLastLevel) ;
06343
06344
06345 for(k=0; k<nn; k++){
06346
06347
06348 for (j = 0; j< (nadjnodnod[k]);j++){
06349
06350
06351 }
06352
06353 priority[k] = W1*(vzdalenosti[k]+1) - W2*(nadjnodnod[k]-1);
06354 statusy[k] = INACTIVE;
06355 priorFronta[k] = 0;
06356
06357
06358 }
06359
06360
06361
06362 oznacenych=0;
06363 statusy[start] = PREACTIVE;
06364 vrcholFronty = 0;
06365 priorFronta[vrcholFronty] = start;
06366
06367 int f=0;
06368
06369
06370 while(vrcholFronty > -1){
06371
06372
06373
06374
06375 for( k=0; k<nn; k++){
06376
06377
06378 }
06379
06380
06381 f=0;
06382 while(f<=vrcholFronty){
06383
06384 f++;
06385 }
06386
06387
06388
06389 maxPriorita = -9999999;
06390
06391 for(k=0; k<=vrcholFronty; k++){
06392
06393
06394 if(priority[priorFronta[k]] > maxPriorita){
06395 maxPriorita = priority[priorFronta[k]];
06396 zpracUzel = priorFronta[k];
06397 poziceUzlu = k;
06398
06399
06400 }
06401
06402
06403 }
06404
06405
06406
06407 priorFronta[poziceUzlu] = priorFronta[vrcholFronty];
06408 vrcholFronty--;
06409 if(statusy[zpracUzel] == PREACTIVE){
06410
06411 for(j = 0; j < nadjnodnod[zpracUzel]; j++){
06412
06413 pomUzel = adjnodnod[zpracUzel][j];
06414 if(pomUzel == zpracUzel)continue;
06415
06416 priority[pomUzel] = priority[pomUzel] + W2;
06417 if(statusy[pomUzel] == INACTIVE){
06418
06419 vrcholFronty++;
06420 priorFronta[vrcholFronty] = pomUzel;
06421 statusy[pomUzel] = PREACTIVE;
06422
06423
06424 }
06425 }
06426
06427 }
06428
06429
06430
06431 oznacenych++;
06432 noveCislovani[zpracUzel] = oznacenych;
06433 statusy[zpracUzel] = POSTACTIVE;
06434
06435 for(j = 0; j < nadjnodnod[zpracUzel]; j++){
06436
06437
06438 pomUzel = adjnodnod[zpracUzel][j];
06439
06440 if(pomUzel == zpracUzel)continue;
06441 if(statusy[pomUzel] == PREACTIVE){
06442
06443
06444
06445 statusy[pomUzel] = ACTIVE;
06446 priority[pomUzel] = priority[pomUzel] + W2;
06447
06448
06449 for(k = 0; k < nadjnodnod[pomUzel]; k++){
06450
06451 pomActive = adjnodnod[pomUzel][k];
06452 if(pomUzel == pomActive)continue;
06453
06454 if(statusy[pomActive] != POSTACTIVE){
06455 priority[pomActive] = priority[pomActive] + W2;
06456 }
06457 if(statusy[pomActive] == INACTIVE){
06458
06459 vrcholFronty++;
06460 priorFronta[vrcholFronty] = pomActive;
06461 statusy[pomActive] = PREACTIVE;
06462
06463
06464 }
06465 }
06466 }
06467
06468 }
06469 }
06470
06471
06472
06473 for( int k=0; k<nn; k++){
06474 ordering[noveCislovani[k]-1] = k;
06475
06476
06477 }
06478
06479
06480
06481
06482
06483
06484
06485
06486
06487 }
06488
06489
06490
06491
06492
06493
06494
06495
06496
06497
06498
06499
06500
06501
06502
06503
06504
06505 void gtopology::searching_hanging_nodes (FILE *out)
06506 {
06507 long i,j,k,nodid,nne,auxnne,ndofn,ndofe;
06508 ivector nod;
06509
06510
06511 for (i=0;i<ne;i++){
06512
06513 nne = give_original_nne (i);
06514 allocv (nne,nod);
06515
06516 give_original_nodes (i,nod);
06517
06518 auxnne=0;
06519
06520 for (j=0;j<nne;j++){
06521
06522 ndofn = gnodes[nod[j]].ndofn;
06523 if (ndofn<0){
06524
06525
06526 auxnne-=ndofn;
06527 }else{
06528 auxnne++;
06529 }
06530 }
06531
06532 if (auxnne>nne){
06533
06534
06535
06536 hangnod=1;
06537
06538
06539
06540 gelements[i].master_nodes = new long [auxnne];
06541 gelements[i].nmne=auxnne;
06542
06543 auxnne=0;
06544 ndofe=0;
06545 for (j=0;j<nne;j++){
06546 nodid=nod[j];
06547 ndofn = gnodes[nodid].ndofn;
06548 if (ndofn>=0){
06549
06550 gelements[i].master_nodes[auxnne]=nodid;
06551 ndofe+=give_ndofn (nodid);
06552 auxnne++;
06553 }
06554 if (ndofn<0){
06555
06556 for (k=0;k<0-ndofn;k++){
06557 gelements[i].master_nodes[auxnne]=gnodes[nodid].mnodes[k];
06558 ndofe+=give_ndofn (gnodes[nodid].mnodes[k]);
06559 auxnne++;
06560 }
06561 }
06562 }
06563
06564
06565 gelements[i].ndofe=ndofe;
06566
06567 gelements[i].cn = new long [ndofe];
06568 }
06569
06570 destrv (nod);
06571 }
06572 }
06573
06574
06575
06576
06577
06578
06579
06580
06581
06582 void gtopology::dof_transfer_hangnod (FILE *out)
06583 {
06584 long i,j,k,l,nmne,ndofn;
06585 ivector nod;
06586
06587
06588 for (i=0;i<ne;i++){
06589 if (gelements[i].master_nodes!=NULL){
06590
06591
06592 gelements[i].cne=1;
06593
06594
06595 nmne = give_nmne (i);
06596
06597 allocv (nmne,nod);
06598
06599 give_master_nodes (i,nod);
06600
06601 l=0;
06602
06603 for (j=0;j<nmne;j++){
06604
06605 ndofn = give_ndofn (nod[j]);
06606
06607
06608 for (k=0;k<ndofn;k++){
06609 gelements[i].cn[l]=give_dof (nod[j],k);
06610 l++;
06611 }
06612 }
06613
06614 destrv (nod);
06615 }
06616 }
06617 }
06618
06619
06620
06621
06622
06623
06624
06625
06626
06627
06628
06629
06630
06631
06632
06633
06634 void gtopology::approx_weights (gtypel et, long nid,vector &w, FILE *out)
06635 {
06636 double xi,eta,zeta;
06637
06638 switch (et){
06639 case isolinear1d:{
06640 if (w.n != 2){
06641 print_err("wrong number of components of the vector of weights", __FILE__, __LINE__, __func__);
06642 }
06643 xi = gnodes[nid].natcoord[0];
06644 bf_lin_1d (w.a,xi);
06645 break;
06646 }
06647 case isoquadratic1d:{
06648 if (w.n != 3){
06649 print_err("wrong number of components of the vector of weights", __FILE__, __LINE__, __func__);
06650 }
06651 xi = gnodes[nid].natcoord[0];
06652 bf_quad_1d (w.a,xi);
06653 break;
06654 }
06655 case trianglelinear:{
06656 if (w.n != 3){
06657 print_err("wrong number of components of the vector of weights", __FILE__, __LINE__, __func__);
06658 }
06659 xi = gnodes[nid].natcoord[0];
06660 eta = gnodes[nid].natcoord[1];
06661 bf_lin_3_2d (w.a,xi,eta);
06662 break;
06663 }
06664 case trianglequadratic:{
06665 if (w.n != 6){
06666 print_err("wrong number of components of the vector of weights", __FILE__, __LINE__, __func__);
06667 }
06668 xi = gnodes[nid].natcoord[0];
06669 eta = gnodes[nid].natcoord[1];
06670 bf_quad_3_2d (w.a,xi,eta);
06671 break;
06672 }
06673 case isolinear2d:{
06674 if (w.n != 4){
06675 print_err("wrong number of components of the vector of weights", __FILE__, __LINE__, __func__);
06676 }
06677 xi = gnodes[nid].natcoord[0];
06678 eta = gnodes[nid].natcoord[1];
06679 bf_lin_4_2d (w.a,xi,eta);
06680 break;
06681 }
06682 case isoquadratic2d:{
06683 if (w.n != 8){
06684 print_err("wrong number of components of the vector of weights", __FILE__, __LINE__, __func__);
06685 }
06686 xi = gnodes[nid].natcoord[0];
06687 eta = gnodes[nid].natcoord[1];
06688 bf_quad_4_2d (w.a,xi,eta);
06689 break;
06690 }
06691 case tetrahedronlinear:{
06692 if (w.n != 4){
06693 print_err("wrong number of components of the vector of weights", __FILE__, __LINE__, __func__);
06694 }
06695 xi = gnodes[nid].natcoord[0];
06696 eta = gnodes[nid].natcoord[1];
06697 zeta = gnodes[nid].natcoord[2];
06698 bf_lin_tet (w.a,xi,eta,zeta);
06699 break;
06700 }
06701 case tetrahedronquadratic:{
06702 if (w.n != 10){
06703 print_err("wrong number of components of the vector of weights", __FILE__, __LINE__, __func__);
06704 }
06705 xi = gnodes[nid].natcoord[0];
06706 eta = gnodes[nid].natcoord[1];
06707 zeta = gnodes[nid].natcoord[2];
06708 bf_quad_tet (w.a,xi,eta,zeta);
06709 break;
06710 }
06711 case isolinear3d:{
06712 if (w.n != 8){
06713 print_err("wrong number of components of the vector of weights", __FILE__, __LINE__, __func__);
06714 }
06715 xi = gnodes[nid].natcoord[0];
06716 eta = gnodes[nid].natcoord[1];
06717 zeta = gnodes[nid].natcoord[2];
06718 bf_lin_hex_3d (w.a,xi,eta,zeta);
06719 break;
06720 }
06721 case isoquadratic3d:{
06722 if (w.n != 20){
06723 print_err("wrong number of components of the vector of weights", __FILE__, __LINE__, __func__);
06724 }
06725 xi = gnodes[nid].natcoord[0];
06726 eta = gnodes[nid].natcoord[1];
06727 zeta = gnodes[nid].natcoord[2];
06728 bf_quad_hex_3d (w.a,xi,eta,zeta);
06729 break;
06730 }
06731 default:{
06732 print_err("unknown type of entity is required", __FILE__, __LINE__, __func__);
06733 }
06734 }
06735 }
06736
06737
06738
06739
06740
06741
06742
06743
06744 void gtopology::hang_nodes_check (FILE *out)
06745 {
06746 long i,j,ndofn,mnid;
06747 double xm,ym,zm,zero;
06748 vector w,x,y,z;
06749
06750 zero=1.0e-2;
06751
06752
06753 for (i=0;i<nn;i++){
06754
06755 ndofn = gnodes[i].give_ndofn ();
06756 if (ndofn<0){
06757
06758
06759 ndofn=0-ndofn;
06760 reallocv (ndofn,w);
06761 reallocv (ndofn,x);
06762 reallocv (ndofn,y);
06763 reallocv (ndofn,z);
06764
06765
06766 approx_weights (gnodes[i].masentity,i,w,out);
06767
06768
06769 for (j=0;j<ndofn;j++){
06770
06771 mnid = gnodes[i].mnodes[j];
06772
06773
06774 x[j] = gnodes[mnid].x;
06775 y[j] = gnodes[mnid].y;
06776 z[j] = gnodes[mnid].z;
06777 }
06778
06779
06780 scprd(x,w,xm);
06781 scprd(y,w,ym);
06782 scprd(z,w,zm);
06783
06784 if (fabs(gnodes[i].x-xm)>zero){
06785 print_err("computed x-coordinate of the hanging node %ld is different from coordinate from input file (%le x %le)",
06786 __FILE__,__LINE__,__func__,i+1,xm,gnodes[i].x);
06787 }
06788 if (fabs(gnodes[i].y-ym)>zero){
06789 print_err("computed y-coordinate of the hanging node %ld is different from coordinate from input file (%le x %le)",
06790 __FILE__,__LINE__,__func__,i+1,ym,gnodes[i].y);
06791 }
06792 if (fabs(gnodes[i].z-zm)>zero){
06793 print_err("computed z-coordinate of the hanging node %ld is different from coordinate from input file (%le x %le)",
06794 __FILE__,__LINE__,__func__,i+1,zm,gnodes[i].z);
06795 }
06796 }
06797 }
06798 }
06799
06800
06801
06802
06803
06804
06805
06806
06807
06808
06809
06810
06811
06812 void gtopology::code_num_mod (long ncomp,long &ndof,FILE *out)
06813 {
06814 long i,j,k,ndofe;
06815 long *cnadd;
06816 ivector cn;
06817 cnadd = new long [ncomp];
06818
06819 for (i=0;i<ncomp;i++){
06820 ndof++;
06821 cnadd[i]=ndof;
06822 }
06823
06824 for (i=0;i<ne;i++){
06825 ndofe = give_ndofe (i)-ncomp;
06826 reallocv (ndofe,cn);
06827 give_code_numbers (i,cn.a);
06828
06829 gelements[i].cne=1;
06830 gelements[i].cn = new long [ndofe+ncomp];
06831 for (j=0;j<ndofe;j++){
06832 gelements[i].cn[j] = cn.a[j];
06833 }
06834 k=0;
06835 for (j=ndofe;j<ndofe+ncomp;j++){
06836 gelements[i].cn[j] = cnadd[k];
06837 k++;
06838 }
06839 }
06840
06841 delete [] cnadd;
06842 }