00001 #include <string.h>
00002 #include <limits.h>
00003 #include <math.h>
00004 #include <float.h>
00005 #include "sequent.h"
00006 #include "global.h"
00007 #include "matrix.h"
00008 #include "vector.h"
00009 #include "gtopology.h"
00010 #include "mathem.h"
00011 #include "element.h"
00012 #include "elemhead.h"
00013 #include "intpoints.h"
00014 #ifndef M_PI
00015 #define M_PI 3.14159265358979323846
00016 #endif
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032 void adjacip (void)
00033 {
00034 time_t bt,et;
00035 bt = time (NULL);
00036
00037 long i,j,k,l,ii,jj,nb,ipp,nip,totnip,nae,nnae,nte,nlmid;
00038 long *adjelel,*adjelelback,*treatedel,*change;
00039 double rlim;
00040
00041
00042 vector coord(3),auxcoord(3);
00043
00044 long adjelel_max = 0;
00045 long adjelelback_max = 0;
00046 long treatedel_max = 0;
00047 long change_max = 0;
00048 long max_nadjip = 0;
00049 long nipp = 0;
00050
00051 adjelel = NULL;
00052 adjelelback = NULL;
00053 treatedel = NULL;
00054 change = NULL;
00055
00056 totnip=Mm->tnip;
00057 Mt->nadjip = new long [totnip];
00058 Mt->adjip = new long* [totnip];
00059 Mt->dist = new double* [totnip];
00060 memset(Mt->nadjip, 0, sizeof(*Mt->nadjip)*totnip);
00061 memset(Mt->adjip, 0, sizeof(*Mt->adjip)*totnip);
00062 memset(Mt->dist, 0, sizeof(*Mt->dist)*totnip);
00063
00064 fprintf(stdout, "\n Adjacent integration points\n");
00065
00066
00067 for (i=0;i<Gtm->ne;i++){
00068 nb = Mt->give_nb (i);
00069 if(i%1000 == 0)
00070 fprintf(stdout, "\n Adjacent integration points on element %ld, max_nadjip=%ld", i+1,max_nadjip);
00071 for (ii=0;ii<nb;ii++){
00072 for (jj=0;jj<nb;jj++){
00073
00074 ipp=Mt->elements[i].ipp[ii][jj];
00075 if (!(Mm->ip[ipp].hmt & 2)){
00076
00077 Mt->nadjip[ipp]=0;
00078 continue;
00079 }
00080 nip = Mt->give_nip (i,ii,jj);
00081
00082 for (j=0;j<nip;j++){
00083
00084 if (!(Mm->ip[ipp].hmt & 2)){
00085
00086 continue;
00087 }
00088
00089 nipp++;
00090
00091 nlmid = Mm->givenonlocid(ipp);
00092
00093 rlim = Mm->nonlocradius (ipp,nlmid);
00094
00095
00096 ipcoord (i,ipp,ii,jj,coord);
00097
00098
00099 Mt->nadjip[ipp]=0;
00100
00101 nae=Gtm->nadjelel[i];
00102 nte=nae;
00103
00104 reallocate(adjelel, nae, adjelel_max);
00105 reallocate(treatedel, nae, treatedel_max);
00106 for (k=0;k<nae;k++){
00107 adjelel[k]=Gtm->adjelel[i][k];
00108 treatedel[k]=adjelel[k];
00109 }
00110
00111
00112
00113
00114 do{
00115 in_dist (ipp,ii,jj,coord,adjelel,nae,rlim);
00116
00117 reallocate(adjelelback, nae, adjelelback_max);
00118
00119 newnadjelel (ipp,adjelelback,adjelel,nae,nnae);
00120 if (nnae>0){
00121
00122 reallocate(adjelel, nnae, adjelel_max);
00123
00124 newadjelel (adjelelback,adjelel,nae,nnae,treatedel,nte);
00125
00126 reallocate(change, nte, change_max);
00127
00128 for (k=0;k<nte;k++){
00129 change[k]=treatedel[k];
00130 }
00131
00132 reallocate(treatedel, nte+nnae, treatedel_max);
00133
00134 for (k=0;k<nte;k++){
00135 treatedel[k]=change[k];
00136 }
00137 l=0;
00138 for (k=nte;k<nte+nnae;k++){
00139 treatedel[k]=adjelel[l]; l++;
00140 }
00141 nte+=nnae;
00142 nae=nnae;
00143 }
00144 else{}
00145 }
00146 while (nnae!=0);
00147
00148
00149
00150
00151
00152 if (max_nadjip < Mt->nadjip[ipp])
00153 max_nadjip = Mt->nadjip[ipp];
00154 Mt->adjip[ipp] = new long [Mt->nadjip[ipp]];
00155 Mt->dist[ipp] = new double [Mt->nadjip[ipp]];
00156
00157
00158 Mt->nadjip[ipp]=0;
00159 nae=Gtm->nadjelel[i];
00160 nte=nae;
00161 reallocate(adjelel, nae, adjelel_max);
00162 reallocate(treatedel, nae, treatedel_max);
00163 for (k=0;k<nae;k++){
00164 adjelel[k]=Gtm->adjelel[i][k];
00165 treatedel[k]=adjelel[k];
00166 }
00167
00168 do{
00169 dist (ipp,ii,jj,max_nadjip,coord,adjelel,nae,rlim);
00170
00171 reallocate(adjelelback, nae, adjelelback_max);
00172
00173 newnadjelel (ipp,adjelelback,adjelel,nae,nnae);
00174 if (nnae>0){
00175
00176 reallocate(adjelel, nnae, adjelel_max);
00177
00178 newadjelel (adjelelback,adjelel,nae,nnae,treatedel,nte);
00179
00180 reallocate(change, nte, change_max);
00181
00182 for (k=0;k<nte;k++){
00183 change[k]=treatedel[k];
00184 }
00185
00186 reallocate(treatedel, nte+nnae, treatedel_max);
00187
00188 for (k=0;k<nte;k++){
00189 treatedel[k]=change[k];
00190 }
00191 l=0;
00192 for (k=nte;k<nte+nnae;k++){
00193 treatedel[k]=adjelel[l]; l++;
00194 }
00195 nte+=nnae;
00196 nae=nnae;
00197 }
00198 }
00199 while (nnae!=0);
00200
00201 ipp++;
00202 }
00203 }
00204 }
00205 }
00206
00207 delete [] adjelelback;
00208 delete [] adjelel;
00209 delete [] treatedel;
00210 delete [] change;
00211
00212
00213
00214
00215
00216
00217
00218
00219
00220
00221
00222
00223
00224
00225 et = time (NULL);
00226 fprintf (Out,"\n\n time of adajcip %ld",et-bt);
00227 fprintf (stdout,"\n\n time of adajcip %ld",et-bt);
00228 }
00229
00230
00231
00232
00233
00234
00235
00236
00237
00238
00239
00240 void save_adjacip()
00241 {
00242 FILE *aux;
00243 char name[FNAMELEN+20];
00244 char emsg[FNAMELEN+100];
00245
00246 switch(Mp->hdbcont.hdbtype)
00247 {
00248 case hdbs_nonloc:
00249 case hdbrs_nonloc:
00250 case hdbrs_single:
00251 case hdbrs_multiple:
00252 case hdbs_single:
00253 case hdbs_multiple:
00254 sprintf(name, "%s.adjacip.bac",Mp->hdbcont.hdbnames);
00255 if (Mp->hdbcont.hdbfmts == text)
00256 aux = fopen(name,"wt");
00257 else
00258 aux = fopen(name,"wb");
00259 if (aux==NULL)
00260 {
00261 sprintf(emsg, "cannot open backup file %s", name);
00262 print_err(emsg, __FILE__, __LINE__, __func__);
00263 abort();
00264 }
00265 if (Mp->hdbcont.hdbfmts == text)
00266 save_adjacip_txt(aux);
00267 else
00268 save_adjacip_bin(aux);
00269 break;
00270 default:
00271 break;
00272 }
00273 }
00274
00275
00276
00277
00278
00279
00280
00281
00282
00283
00284
00285
00286 long restore_adjacip()
00287 {
00288 FILE *aux;
00289 char name[FNAMELEN+20];
00290 char emsg[FNAMELEN+100];
00291
00292 switch(Mp->hdbcont.hdbtype)
00293 {
00294 case hdbr_nonloc:
00295 case hdbrs_nonloc:
00296 case hdbrs_single:
00297 case hdbrs_multiple:
00298 case hdbr_single:
00299 case hdbr_multiple:
00300 sprintf(name, "%s.adjacip.bac",Mp->hdbcont.hdbnamer);
00301 if (Mp->hdbcont.hdbfmtr == text)
00302 aux = fopen(name,"rt");
00303 else
00304 aux = fopen(name,"rb");
00305 if (aux==NULL)
00306 {
00307 sprintf(emsg, "cannot open backup file %s", name);
00308 print_err(emsg, __FILE__, __LINE__, __func__);
00309 abort();
00310 }
00311 if (Mp->hdbcont.hdbfmtr == text)
00312 restore_adjacip_txt(aux);
00313 else
00314 restore_adjacip_bin(aux);
00315 break;
00316 default:
00317 return 0;
00318 }
00319 return 1;
00320 }
00321
00322
00323
00324
00325
00326
00327
00328
00329
00330
00331
00332
00333
00334 void save_adjacip_txt(FILE *aux)
00335 {
00336 long i,j;
00337 if (Mespr)
00338 fprintf(stdout, "\n Saving adjacent integration points to backup file");
00339 for (i=0;i<Mm->tnip;i++)
00340 {
00341 fprintf (aux,"%ld %ld",i,Mt->nadjip[i]);
00342 for (j=0;j<Mt->nadjip[i];j++)
00343 fprintf (aux," %ld %le",Mt->adjip[i][j],Mt->dist[i][j]);
00344 fprintf(aux, "\n");
00345 }
00346 }
00347
00348
00349
00350
00351
00352
00353
00354
00355
00356
00357
00358
00359
00360 void save_adjacip_bin(FILE *aux)
00361 {
00362 long i,j;
00363 if (Mespr)
00364 fprintf(stdout, "\n Saving adjacent integration points to backup file");
00365 for (i=0;i<Mm->tnip;i++)
00366 {
00367 fwrite(&i, sizeof(i), 1, aux);
00368 fwrite(Mt->nadjip+i, sizeof(*Mt->nadjip), 1, aux);
00369 for (j=0;j<Mt->nadjip[i];j++)
00370 {
00371 fwrite(Mt->adjip[i]+j, sizeof(*Mt->adjip[i]), 1, aux);
00372 fwrite(Mt->dist[i]+j, sizeof(*Mt->dist[i]), 1, aux);
00373 }
00374 }
00375 }
00376
00377
00378
00379
00380
00381
00382
00383
00384
00385
00386
00387
00388
00389 void restore_adjacip_txt(FILE *aux)
00390 {
00391 long i,j,ipp;
00392 Mt->nadjip = new long[Mm->tnip];
00393 memset(Mt->nadjip, 0, sizeof(*Mt->nadjip)*Mm->tnip);
00394 Mt->adjip = new long*[Mm->tnip];
00395 memset(Mt->adjip, 0, sizeof(*Mt->adjip)*Mm->tnip);
00396 Mt->dist = new double*[Mm->tnip];
00397 memset(Mt->dist, 0, sizeof(*Mt->dist)*Mm->tnip);
00398 if (Mespr)
00399 fprintf(stdout, "\n Restoring adjacent integration points from backup file");
00400 for (i=0;i<Mm->tnip;i++)
00401 {
00402 fscanf (aux, "%ld%ld", &ipp, Mt->nadjip+i);
00403 if ((ipp < 0) || (ipp >= Mm->tnip))
00404 {
00405 print_err("invalid integration point number", __FILE__, __LINE__, __func__);
00406 abort();
00407 }
00408 Mt->adjip[i] = new long[Mt->nadjip[i]];
00409 memset(Mt->adjip[i], 0, sizeof(*Mt->adjip[i])*Mt->nadjip[i]);
00410 Mt->dist[i] = new double[Mt->nadjip[i]];
00411 memset(Mt->dist[i], 0, sizeof(*Mt->dist[i])*Mt->nadjip[i]);
00412 for (j=0; j<Mt->nadjip[i]; j++)
00413 fscanf (aux, "%ld%le", Mt->adjip[i]+j, Mt->dist[i]+j);
00414 }
00415 }
00416
00417
00418
00419
00420
00421
00422
00423
00424
00425
00426
00427
00428
00429 void restore_adjacip_bin(FILE *aux)
00430 {
00431 long i,j,ipp;
00432 Mt->nadjip = new long[Mm->tnip];
00433 memset(Mt->nadjip, 0, sizeof(*Mt->nadjip)*Mm->tnip);
00434 Mt->adjip = new long*[Mm->tnip];
00435 memset(Mt->adjip, 0, sizeof(*Mt->adjip)*Mm->tnip);
00436 Mt->dist = new double*[Mm->tnip];
00437 memset(Mt->dist, 0, sizeof(*Mt->dist)*Mm->tnip);
00438 if (Mespr)
00439 fprintf(stdout, "\n Restoring adjacent integration points from backup file");
00440 for (i=0;i<Mm->tnip;i++)
00441 {
00442 fread(&ipp, sizeof(i), 1, aux);
00443 fread(Mt->nadjip+i, sizeof(*Mt->nadjip), 1, aux);
00444 if ((ipp < 0) || (ipp >= Mm->tnip))
00445 {
00446 print_err("invalid integration point number", __FILE__, __LINE__, __func__);
00447 abort();
00448 }
00449 Mt->adjip[i] = new long[Mt->nadjip[i]];
00450 memset(Mt->adjip[i], 0, sizeof(*Mt->adjip[i])*Mt->nadjip[i]);
00451 Mt->dist[i] = new double[Mt->nadjip[i]];
00452 memset(Mt->dist[i], 0, sizeof(*Mt->dist[i])*Mt->nadjip[i]);
00453 for (j=0;j<Mt->nadjip[i];j++)
00454 {
00455 fread(Mt->adjip[i]+j, sizeof(*Mt->adjip[i]), 1, aux);
00456 fread(Mt->dist[i]+j, sizeof(*Mt->dist[i]), 1, aux);
00457 }
00458 }
00459 }
00460
00461
00462
00463
00464
00465
00466
00467
00468
00469
00470 long give_max_adjacip()
00471 {
00472 long i,tnip,max_tnip,ipp,ret,nlmid;
00473 long max_1d, max_2d, max_3d;
00474 double vol, area, length, min_vol, min_area, min_length;
00475 double r, max_r;
00476
00477 ret = 0;
00478 max_tnip = 0;
00479 max_r = 0.0;
00480 min_vol = min_area = min_length = DBL_MAX;
00481 for (i=0; i<Mt->ne; i++)
00482 {
00483 ipp=Mt->elements[i].ipp[0][0];
00484
00485 if (!(Mm->ip[ipp].hmt & 2))
00486
00487 continue;
00488
00489 nlmid = Mm->givenonlocid(ipp);
00490
00491 r = Mm->nonlocradius (ipp,nlmid);
00492 if (r > max_r)
00493 max_r = r;
00494
00495 tnip = Mt->give_tnip(i);
00496 if (tnip > max_tnip)
00497 max_tnip = tnip;
00498
00499 switch (Mt->give_dimension(i))
00500 {
00501 case 1:
00502 length = fabs(Mt->give_length(i));
00503 if (length < min_length)
00504 min_length = length;
00505 break;
00506 case 2:
00507 area = fabs(Mt->give_area(i));
00508 if (area < min_area)
00509 min_area = area;
00510 break;
00511 case 3:
00512 vol = fabs(Mt->give_volume(i));
00513 if (vol < min_vol)
00514 min_vol = vol;
00515 break;
00516 default:
00517 print_err("unknown dimension of element", __FILE__, __LINE__, __func__);
00518 }
00519 }
00520 if (min_length == 0.0)
00521 print_err("minimum length of element is zero ", __FILE__, __LINE__, __func__);
00522 if (min_area == 0.0)
00523 print_err("minimum area of element is zero ", __FILE__, __LINE__, __func__);
00524 if (min_vol == 0.0)
00525 print_err("minimum volume of element is zero ", __FILE__, __LINE__, __func__);
00526
00527
00528 max_1d = 15*(long(2.0*max_r/min_length));
00529
00530 max_2d = 15*(long(M_PI*max_r*max_r/min_area));
00531
00532 max_3d = 15*(long(4.0/3.0*M_PI*max_r*max_r*max_r/min_vol));
00533 if (max_1d > max_2d)
00534 ret = max_1d;
00535 else
00536 ret = max_2d;
00537 if (max_3d > ret)
00538 ret = max_3d;
00539
00540 return ret;
00541 }
00542
00543
00544
00545
00546
00547
00548
00549
00550
00551
00552
00553
00554
00555
00556
00557 void reallocate(long *&array, long req_size, long &capacity)
00558 {
00559 if (req_size == 0)
00560 return;
00561 if (req_size > capacity)
00562 {
00563
00564 capacity = 2*req_size;
00565 if (array)
00566 delete [] array;
00567 array = new long [capacity];
00568 }
00569 }
00570
00571
00572
00573
00574
00575
00576
00577
00578
00579
00580
00581
00582
00583
00584
00585
00586
00587 void newnadjelel (long ipp,long *adjelelback,long *adjelel,long &nae,long &nnae)
00588 {
00589 long i,nee;
00590
00591
00592 nee=0;
00593 for (i=0;i<nae;i++){
00594 if (adjelel[i]<0) nee++;
00595 }
00596
00597 if (nae!=nee){
00598 nnae=0;
00599 for (i=0;i<nae;i++){
00600 adjelelback[i]=adjelel[i];
00601 if (adjelel[i]>-1){
00602 nnae+=Gtm->nadjelel[adjelel[i]];
00603 }
00604 }
00605 }
00606 else{
00607 nnae=0;
00608 }
00609 }
00610
00611
00612
00613
00614
00615
00616
00617
00618
00619
00620
00621
00622
00623
00624
00625
00626
00627
00628 void newadjelel (long *adjelelback,long *adjelel,long &nae,long &nnae,long *treatedel,long &nte)
00629 {
00630 long i,j,k,l,prev,min;
00631
00632 j=0;
00633 for (i=0;i<nae;i++){
00634 if (adjelelback[i]>-1){
00635 for (k=0;k<Gtm->nadjelel[adjelelback[i]];k++){
00636 adjelel[j]=Gtm->adjelel[adjelelback[i]][k];
00637 j++;
00638 }
00639 }
00640 }
00641
00642
00643 prev=-1;
00644 for (i=0;i<nnae;i++){
00645 min=LONG_MAX;
00646 for (j=i;j<nnae;j++){
00647 if (min>adjelel[j]){
00648 min=adjelel[j]; k=j;
00649 }
00650 }
00651 if (min==prev){
00652 nnae--; i--;
00653 adjelel[k]=adjelel[nnae];
00654 }
00655 else{
00656 l=adjelel[i];
00657 adjelel[i]=min;
00658 adjelel[k]=l;
00659 prev=min;
00660 }
00661 }
00662
00663 for (i=0;i<nnae;i++){
00664 k=adjelel[i];
00665 for (j=0;j<nte;j++){
00666 if (treatedel[j]==k){
00667 nnae--;
00668 adjelel[i]=adjelel[nnae];
00669 i--;
00670 break;
00671 }
00672 }
00673 }
00674
00675 }
00676
00677
00678
00679
00680
00681
00682
00683
00684
00685
00686
00687
00688
00689
00690
00691
00692
00693
00694
00695
00696 void in_dist (long ipp,long ri,long ci,vector &coord,long *adjelel,long &nae,double rlim)
00697 {
00698 long i,j,lj,uj,ii,eid,nip;
00699 double r;
00700 vector auxcoord(3);
00701
00702 for (i=0;i<nae;i++){
00703 eid = adjelel[i];
00704 nip = Mt->give_nip (eid,ri,ci);
00705 lj = Mt->elements[eid].ipp[ri][ci];
00706 uj = lj+nip;
00707
00708 ii=0;
00709 for (j=lj;j<uj;j++){
00710 if (Mm->ip[j].hmt & 2){
00711
00712 ipcoord (eid,j,ri,ci,auxcoord);
00713 r = length (coord,auxcoord);
00714
00715 if (r<rlim){
00716 Mt->nadjip[ipp]++;
00717 ii++;
00718 }
00719 }
00720 }
00721
00722 if (ii==0) adjelel[i]=-1-adjelel[i];
00723
00724 }
00725 }
00726
00727
00728
00729
00730
00731
00732
00733
00734
00735
00736
00737
00738
00739
00740
00741
00742
00743
00744
00745
00746
00747 void dist (long ipp,long ri,long ci,long max_nadjip,vector &coord,long *adjelel,long &nae,double rlim)
00748 {
00749 long i,j,lj,uj,ii,eid,nip;
00750 double r;
00751 vector auxcoord(3);
00752
00753 for (i=0;i<nae;i++){
00754 eid = adjelel[i];
00755 nip = Mt->give_nip (eid,ri,ci);
00756 lj = Mt->elements[eid].ipp[ri][ci];
00757 uj = lj+nip;
00758
00759 ii=0;
00760 for (j=lj;j<uj;j++){
00761 if (Mm->ip[j].hmt & 2){
00762
00763
00764 ipcoord (eid,j,ri,ci,auxcoord);
00765 r = length (coord,auxcoord);
00766
00767 if (r<rlim){
00768 Mt->adjip[ipp][Mt->nadjip[ipp]]=j;
00769 Mt->dist[ipp][Mt->nadjip[ipp]]=r;
00770 Mt->nadjip[ipp]++;
00771 ii++;
00772 }
00773 }
00774 }
00775 if (ii==0) adjelel[i]=-1-adjelel[i];
00776 }
00777 }
00778
00779
00780
00781
00782
00783
00784
00785
00786
00787
00788
00789
00790
00791
00792
00793
00794
00795
00796 void ipcoord (long eid,long ipp,long ri,long ci,vector &ipcoord)
00797 {
00798 switch (Mt->give_elem_type(eid)){
00799 case bar3d:
00800 Bar3d->ipcoord (eid,ipp,ri,ci,ipcoord);
00801 break;
00802 case axisymmlt:
00803 Asymlt->ipcoord(eid, ipp, ri, ci, ipcoord);
00804 break;
00805 case axisymmlq:
00806 Asymlq->ipcoord(eid, ipp, ri, ci, ipcoord);
00807 break;
00808 case axisymmqq:
00809 Asymqq->ipcoord(eid, ipp, ri, ci, ipcoord);
00810 break;
00811 case planeelementlt:{
00812 Pelt->ipcoord (eid,ipp,ri,ci,ipcoord);
00813 break;
00814 }
00815 case planeelementqt:{
00816 Peqt->ipcoord (eid,ipp,ri,ci,ipcoord);
00817 break;
00818 }
00819 case planeelementrotlt:{
00820 Perlt->ipcoord (eid,ipp,ri,ci,ipcoord);
00821 break;
00822 }
00823 case planeelementlq:{
00824 Pelq->ipcoord (eid,ipp,ri,ci,ipcoord);
00825 break;
00826 }
00827 case planeelementqq:{
00828 Peqq->ipcoord (eid,ipp,ri,ci,ipcoord);
00829 break;
00830 }
00831 case planeelementrotlq:{
00832 Perlq->ipcoord (eid,ipp,ri,ci,ipcoord);
00833 break;
00834 }
00835
00836 case lineartet:{
00837 Ltet->ipcoord (eid,ipp,ri,ci,ipcoord);
00838 break;
00839 }
00840 case quadrtet:{
00841 Qtet->ipcoord (eid,ipp,ri,ci,ipcoord);
00842 break;
00843 }
00844 case linearhex:{
00845 Lhex->ipcoord (eid,ipp,ri,ci,ipcoord);
00846 break;
00847 }
00848 case quadrhex:{
00849 Qhex->ipcoord (eid,ipp,ri,ci,ipcoord);
00850 break;
00851 }
00852 default:{
00853 print_err("unknown element type is required", __FILE__, __LINE__, __func__);
00854 }
00855 }
00856 }