00001 #include "loadcaset.h"
00002 #include "globalt.h"
00003 #include "aliast.h"
00004 #include "globmatt.h"
00005 #include "genfile.h"
00006 #include "elemswitcht.h"
00007 #include <math.h>
00008 #include <limits.h>
00009
00010 loadcaset::loadcaset (void)
00011 {
00012
00013 neb = 0;
00014
00015
00016 nqs = 0;
00017
00018 nnqs = 0;
00019
00020 lnqs = NULL;
00021
00022 neqs = 0;
00023
00024 leqs = NULL;
00025
00026 neins=0;
00027
00028 en = NULL;
00029
00030 sip = 0;
00031
00032 lnc = NULL;
00033
00034 nc = NULL;
00035
00036 ne = 0;
00037
00038 sour = NULL;
00039
00040
00041 nnv=0;
00042
00043 nodval=NULL;
00044
00045 ncc=0;
00046 ncc2=0;
00047
00048 climcond=NULL;
00049 climcond2=NULL;
00050
00051
00052
00053
00054 npv = 0;
00055
00056 pv = NULL;
00057
00058 elemload = NULL;
00059 climcond = NULL;
00060 climcond2 = NULL;
00061 tabf = NULL;
00062
00063
00064 masterval = 0.0;
00065 mastergrad = NULL;
00066 ndim = 0;
00067 }
00068
00069 loadcaset::~loadcaset (void)
00070 {
00071 long i;
00072
00073
00074 for (i=0;i<nnqs;i++){
00075 delete [] lnqs[i];
00076 }
00077 delete [] lnqs;
00078
00079
00080 for (i=0;i<neqs;i++){
00081 delete [] leqs[i];
00082 }
00083 delete [] leqs;
00084
00085
00086 delete en;
00087
00088
00089 if(lnc != NULL){
00090 for (i=0;i<nqs;i++){
00091 delete [] lnc[i];
00092 }
00093 delete [] lnc;
00094 }
00095
00096
00097 delete [] nc;
00098
00099
00100 delete [] sour;
00101
00102 delete [] nodval;
00103
00104
00105 delete [] pv;
00106
00107 delete [] tabf;
00108 delete [] elemload;
00109 delete [] climcond;
00110 delete [] climcond2;
00111
00112 delete [] mastergrad;
00113 }
00114
00115
00116
00117
00118
00119
00120
00121
00122
00123 void loadcaset::read (XFILE *in,long lcid)
00124 {
00125 long i,j,k;
00126 XFILE *locin;
00127 char fname[1020],*path,*name,*suffix;
00128
00129
00130
00131
00132
00133 xfscanf (in,"%ld",&npv);
00134 if (Mesprt==1)
00135 fprintf (stdout,"\n the number of prescribed values %ld",npv);
00136 if (npv<0)
00137 print_err("negative number of prescribed values",__FILE__,__LINE__,__func__);
00138
00139 pv = new pvalt [npv];
00140 for (i=0;i<npv;i++){
00141 pv[i].read(in);
00142 }
00143
00144
00145
00146
00147
00148 xfscanf (in,"%ld",&nqs);
00149 if (Mesprt==1)
00150 fprintf (stdout,"\n number of prescribed sources %ld",nqs);
00151 if (nqs<0)
00152 print_err("negative number of quantity sources",__FILE__,__LINE__,__func__);
00153 if (nqs>0){
00154
00155 sour = new sourcet [nqs];
00156 for (i=0;i<nqs;i++){
00157 xfscanf (in,"%ld",&j);
00158 if (j<1 || j>nqs)
00159 print_err("wrong number of quantity source",__FILE__,__LINE__,__func__);
00160 sour[j-1].read (in);
00161 }
00162
00163
00164 xfscanf (in,"%ld",&nnqs);
00165 if (nnqs<0)
00166 print_err("negative number of nodes with defined source",__FILE__,__LINE__,__func__);
00167
00168 lnqs = new long* [nnqs];
00169 for (i=0;i<nnqs;i++){
00170 lnqs[i] = new long [2];
00171 xfscanf (in,"%ld %ld",&lnqs[i][0],&lnqs[i][1]);
00172 lnqs[i][0]--;
00173 lnqs[i][1]--;
00174 }
00175
00176
00177 xfscanf (in,"%ld",&neqs);
00178 if (neqs<0)
00179 print_err("negative number of elements with defined source",__FILE__,__LINE__,__func__);
00180
00181 if (Mesprt==1)
00182 fprintf (stdout,"\n number of defined sources %ld",neqs);
00183
00184 leqs = new long* [neqs];
00185 for (i=0;i<neqs;i++){
00186 leqs[i] = new long [2];
00187 xfscanf (in,"%ld %ld",&leqs[i][0],&leqs[i][1]);
00188 leqs[i][0]--;
00189 leqs[i][1]--;
00190 }
00191 }
00192
00193
00194
00195
00196
00197 xfscanf (in,"%ld",&neb);
00198 if (neb<0)
00199 print_err("negative number of elements with boundary conditions",__FILE__,__LINE__,__func__);
00200
00201 elemload = new loadelt [neb];
00202 if (Mesprt==1)
00203 fprintf (stdout,"\n number of elements with boundary edges %ld",neb);
00204 for (i=0;i<neb;i++){
00205 elemload[i].read (in,lcid);
00206 }
00207
00208
00209
00210
00211
00212 xfscanf (in,"%ld",&nnv);
00213 if (nnv<0)
00214 print_err("negative number of objects with nodal values",__FILE__,__LINE__,__func__);
00215
00216 if (Mesprt==1)
00217 fprintf (stdout,"\n number of objects storing nodal values %ld",nnv);
00218
00219 nodval = new bnodvalt [nnv];
00220 for (i=0;i<nnv;i++){
00221 nodval[i].read (in);
00222 }
00223
00224
00225
00226
00227
00228 xfscanf (in,"%ld",&ncc);
00229 if (ncc<0)
00230 print_err("negative number of objects with nodal values",__FILE__,__LINE__,__func__);
00231
00232 if (Mesprt==1)
00233 fprintf (stdout,"\n number of objects of climatic conditions %ld",ncc);
00234
00235 climcond = new climatcond [ncc];
00236 if (ncc>0){
00237
00238 xfscanf (in,"%ld",&trcc);
00239
00240 if (trcc==1){
00241
00242 for (i=0;i<ncc;i++){
00243 xfscanf (in,"%ld",&j);
00244 if (j<1 || j>ncc){
00245 print_err("index of climate condition record is out of range",__FILE__,__LINE__,__func__);
00246 }
00247 climcond[j-1].read(in);
00248 }
00249 }
00250 if (trcc==2){
00251
00252
00253
00254 filename_decomposition(in->fname,path,name,suffix);
00255 sprintf(fname, "%s%s.txt", path, name);
00256
00257 locin = xfopen(fname,"r");
00258 delete [] path; delete [] name; delete [] suffix;
00259
00260 locin->warning = in->warning;
00261 locin->kwdmode = in->kwdmode;
00262 locin->ignorecase = in->ignorecase;
00263
00264 for (i=0;i<ncc;i++){
00265 xfscanf (locin,"%ld",&j);
00266 if (j<1 || j>ncc){
00267 print_err("index of climate condition record is out of range",__FILE__,__LINE__,__func__);
00268 }
00269 climcond[j-1].read (locin);
00270 }
00271
00272 xfclose (locin);
00273 }
00274 }
00275
00276
00277 xfscanf (in,"%ld",&ncc2);
00278 if (ncc2<0)
00279 print_err("negative number of objects with nodal values",__FILE__,__LINE__,__func__);
00280
00281 if (Mesprt==1)
00282 fprintf (stdout,"\n number of objects of climatic conditions type 2 %ld",ncc2);
00283
00284 climcond2 = new climatcond2 [ncc2];
00285 if (ncc2>0){
00286
00287 xfscanf (in,"%ld",&trcc2);
00288
00289 if (trcc2==1){
00290
00291 for (i=0;i<ncc2;i++){
00292 xfscanf (in,"%ld",&j);
00293 if (j<1 || j>ncc2){
00294 print_err("index of climate condition record is out of range",__FILE__,__LINE__,__func__);
00295 }
00296 climcond2[j-1].read (in);
00297 }
00298 }
00299 if (trcc2==2){
00300
00301
00302
00303 filename_decomposition(in->fname,path,name,suffix);
00304 sprintf(fname, "%s%s.txt", path, name);
00305
00306 locin = xfopen(fname,"r");
00307 delete [] path; delete [] name; delete [] suffix;
00308
00309 locin->warning = in->warning;
00310 locin->kwdmode = in->kwdmode;
00311 locin->ignorecase = in->ignorecase;
00312
00313 for (i=0;i<ncc2;i++){
00314 xfscanf (locin,"%ld",&j);
00315 if (j<1 || j>ncc2){
00316 print_err("index of climate condition record is out of range",__FILE__,__LINE__,__func__);
00317 }
00318 climcond2[j-1].read (locin);
00319 }
00320
00321 xfclose (locin);
00322 }
00323 }
00324
00325
00326
00327
00328
00329
00330
00331
00332
00333
00334 xfscanf (in,"%ld",&nmultfunc);
00335
00336 if (Mesprt==1)
00337 fprintf (stdout,"\n number of objects of time functions %ld",nmultfunc);
00338
00339 tabf = new tablefunct [nmultfunc];
00340
00341 for(i=0;i<nmultfunc;i++){
00342
00343 xfscanf (in,"%ld",&k);
00344
00345 if ((k-1) != i){
00346
00347 fprintf (stderr,"\n\n wrong number of prescribed time function loadcaset::read (file %s, line %d).\n",__FILE__,__LINE__);
00348 exit(0);
00349 }
00350
00351 xfscanf(in,"%d",(int*)&tfunc);
00352 switch (tfunc){
00353
00354 case tab:{
00355 tabf[i].read(in);
00356 break;
00357 }
00358 default:{
00359 fprintf (stderr,"\n\n unknown function is required in function pvalt::read (file %s, line %d).\n",__FILE__,__LINE__);
00360 }
00361 }
00362 }
00363
00364
00365
00366
00367 if(Tp->homogt == 1){
00368 ndim = Tp->gdim;
00369
00370 mastergrad = new double [ndim];
00371
00372 xfscanf (in,"%lf",&masterval);
00373 for (i=0;i<ndim;i++)
00374 xfscanf (in,"%lf",&mastergrad[i]);
00375 }
00376 }
00377
00378
00379
00380
00381
00382
00383
00384
00385
00386
00387 void loadcaset::print (FILE *out,long lcid)
00388 {
00389 long i;
00390
00391
00392
00393
00394 fprintf(out,"\n## dirichlet b.c.:");
00395 fprintf (out,"\n%ld\n\n",npv);
00396 for (i=0;i<npv;i++){
00397 pv[i].print(out);
00398 }
00399
00400
00401
00402
00403
00404 fprintf(out,"\n## sources:");
00405 fprintf (out,"\n%ld\n\n",nqs);
00406 if (nqs>0){
00407
00408 for (i=0;i<nqs;i++){
00409 fprintf (out,"\n %ld",i+1);
00410 sour[i].print (out);
00411 }
00412
00413
00414 fprintf (out,"\n\n %ld",nnqs);
00415 for (i=0;i<nnqs;i++){
00416 fprintf (out,"\n %ld %ld",lnqs[i][0]+1,lnqs[i][1]+1);
00417 }
00418
00419
00420 fprintf (out,"\n\n %ld",neqs);
00421 for (i=0;i<neqs;i++){
00422 fprintf (out,"\n %ld %ld",leqs[i][0]+1,leqs[i][1]+1);
00423 }
00424 }
00425
00426
00427
00428
00429 fprintf(out,"\n## loaded elements:");
00430 fprintf (out,"\n%ld\n\n",neb);
00431 for (i=0;i<neb;i++){
00432 elemload[i].print (out,lcid);
00433 }
00434
00435
00436
00437
00438 fprintf(out,"\n## nodal values:");
00439 fprintf (out,"\n%ld\n\n",nnv);
00440 for (i=0;i<nnv;i++){
00441 nodval[i].print (out);
00442 }
00443
00444
00445
00446
00447 fprintf(out,"\n## climatic conditions:");
00448 fprintf (out,"\n%ld\n",ncc);
00449 if (ncc>0){
00450 fprintf (out,"%ld\n",trcc);
00451 for (i=0;i<ncc;i++){
00452 fprintf (out,"\n\n %ld",i+1);
00453 climcond[i].print (out);
00454 }
00455 }
00456
00457 fprintf(out,"\n## climatic conditions2:");
00458 fprintf (out,"\n%ld\n",ncc2);
00459 if (ncc>0){
00460 fprintf (out,"%ld\n",trcc);
00461 for (i=0;i<ncc;i++){
00462 fprintf (out,"\n\n %ld",i+1);
00463 climcond2[i].print (out);
00464 }
00465 }
00466
00467
00468
00469
00470
00471 fprintf(out,"\n## time functions:");
00472 fprintf (out,"\n%ld\n\n",nmultfunc);
00473
00474 for(i=0;i<nmultfunc;i++){
00475
00476 fprintf (out,"\n %ld",i+1);
00477
00478 fprintf(out,"\n %d ",(int)tfunc);
00479 switch (tfunc){
00480
00481 case tab:{
00482 tabf[i].print(out);
00483 break;
00484 }
00485 default:{
00486 }
00487 }
00488 }
00489
00490
00491
00492
00493 fprintf(out,"\n## homogenization:\n");
00494 if(Tp->homogt == 1){
00495 fprintf (out,"\n %lf",masterval);
00496 for (i=0;i<ndim;i++)
00497 fprintf (out," %lf",mastergrad[i]);
00498 fprintf (out,"\n");
00499 }
00500
00501 }
00502
00503
00504
00505
00506
00507
00508
00509
00510
00511
00512
00513
00514
00515
00516
00517
00518
00519
00520
00521
00522
00523 void loadcaset::assemble (long lcid,double *rhs)
00524 {
00525 long i,ndofe,eid,*cn;
00526 vector r,f,lv,nodval;
00527 matrix km;
00528
00529
00530
00531
00532 source_contrib (lcid,rhs);
00533
00534
00535
00536
00537
00538 for (i=0;i<neb;i++){
00539
00540 eid=elemload[i].eid;
00541
00542 if(eid < 0){
00543 print_err("element number is less than 0", __FILE__, __LINE__, __func__);
00544 abort();
00545 }
00546
00547 if (Gtt->leso[eid]==1){
00548
00549
00550
00551 ndofe=Tt->give_ndofe (eid);
00552 reallocv (ndofe,lv);
00553 fillv (0.0,lv);
00554
00555 if (Tp->tprob == growing_np_problem || Tp->tprob == growing_np_problem_nonlin){
00556
00557 }
00558
00559
00560
00561
00562 elem_neumann_vector (lv,lcid,eid,i);
00563
00564
00565
00566
00567 elem_newton_vector (lv,lcid,eid,i);
00568
00569
00570
00571 cn = new long [ndofe];
00572 Tt->give_code_numbers (eid,cn);
00573
00574
00575 locglob (rhs,lv.a,cn,ndofe);
00576
00577 delete [] cn;
00578 }
00579 }
00580
00581
00582
00583
00584
00585
00586
00587
00588
00589
00590
00591
00592
00593
00594
00595
00596
00597
00598
00599
00600
00601
00602
00603
00604
00605
00606
00607
00608
00609 }
00610
00611
00612
00613
00614
00615
00616
00617
00618
00619
00620
00621
00622
00623
00624
00625
00626 void loadcaset::assemble_flux (long lcid,double *rhs,long n)
00627 {
00628 long i,ndofe,eid,*cn;
00629 vector r,f,lv,nodval;
00630 matrix km;
00631
00632
00633 nullv (rhs,n);
00634
00635
00636
00637
00638 source_contrib (lcid,rhs);
00639
00640
00641
00642
00643 for (i=0;i<neb;i++){
00644
00645 eid=elemload[i].eid;
00646
00647 if(eid < 0){
00648 print_err("element number is less than 0", __FILE__, __LINE__, __func__);
00649 abort();
00650 }
00651
00652 if (Gtt->leso[eid]==1){
00653
00654
00655
00656 ndofe=Tt->give_ndofe (eid);
00657 allocv (ndofe,lv);
00658 fillv (0.0,lv);
00659
00660 if (Tp->tprob == growing_np_problem || Tp->tprob == growing_np_problem_nonlin){
00661
00662 }
00663
00664
00665
00666
00667 elem_neumann_vector (lv,lcid,eid,i);
00668
00669
00670
00671
00672 elem_transmission_flux (lv,lcid,eid,i);
00673
00674 cn = new long [ndofe];
00675 Gtt->give_code_numbers (eid,cn);
00676
00677
00678 locglob (rhs,lv.a,cn,ndofe);
00679
00680 destrv (lv);
00681 delete [] cn;
00682 }
00683 }
00684 }
00685
00686
00687 double loadcaset::give_fact (double t,long i)
00688 {
00689 return tabf[i-1].getval(t);
00690 }
00691
00692
00693
00694
00695
00696
00697
00698 void loadcaset::elemsource ()
00699 {
00700 long i,j;
00701 long *aux;
00702
00703
00704 aux = new long [nnqs];
00705 for (i=0;i<nnqs;i++){
00706 aux[i]=lnqs[i][0];
00707 }
00708
00709 en = new elemnode;
00710
00711 en->selnode (nnqs,aux);
00712
00713 delete [] aux;
00714
00715
00716 en->elemnodes (Gtt);
00717
00718
00719 neins = en->nie;
00720
00721 for (i=0;i<neins;i++){
00722 j=en->lse[i];
00723 Tt->elements[j].source=1;
00724 }
00725
00726 for (i=0;i<neqs;i++){
00727 j=leqs[i][0];
00728 Tt->elements[j].source=1;
00729 }
00730
00731
00732
00733
00734
00735
00736
00737
00738
00739
00740
00741
00742
00743
00744
00745
00746
00747
00748
00749
00750
00751
00752
00753
00754
00755
00756
00757
00758
00759
00760
00761
00762
00763
00764
00765
00766
00767
00768
00769
00770
00771
00772
00773
00774
00775
00776
00777
00778
00779
00780
00781
00782
00783
00784
00785
00786
00787
00788
00789
00790
00791
00792
00793
00794
00795
00796
00797
00798
00799
00800
00801
00802
00803
00804
00805
00806
00807
00808
00809
00810
00811
00812
00813
00814
00815
00816
00817
00818
00819
00820
00821
00822
00823
00824
00825
00826
00827
00828
00829
00830
00831
00832
00833
00834
00835
00836
00837
00838
00839
00840
00841
00842
00843
00844
00845
00846
00847
00848
00849
00850
00851
00852
00853
00854
00855
00856
00857
00858
00859
00860
00861
00862
00863
00864
00865
00866
00867
00868
00869
00870
00871
00872
00873
00874
00875
00876
00877
00878
00879
00880
00881
00882
00883
00884
00885
00886
00887
00888
00889
00890
00891
00892
00893
00894
00895
00896
00897
00898
00899
00900
00901
00902
00903 }
00904
00905
00906
00907
00908
00909
00910
00911
00912
00913 void loadcaset::sourcenodalvalues (long idse,vector &nodval)
00914 {
00915 long i,j,k,m,nne;
00916
00917
00918 nne = nodval.n;
00919
00920
00921 for (i=0;i<nne;i++){
00922
00923 j=en->elnod[idse][i];
00924 if (j>-1){
00925
00926 k=lnqs[j][1];
00927
00928 m=en->lse[idse];
00929 nodval[i]=sour[k].giveval (m);
00930 }
00931 else{
00932 nodval[i]=0.0;
00933 }
00934 }
00935
00936 }
00937
00938
00939
00940
00941
00942
00943
00944
00945
00946 void loadcaset::source_contrib (long lcid,double *rhs)
00947 {
00948 long i,j,k,nne,ndofe,*cn;
00949 double s;
00950 vector lv,nodval;
00951
00952
00953 if (nnqs>0){
00954 for (i=0;i<neins;i++){
00955
00956
00957
00958 j=en->lse[i];
00959 if (Gtt->leso[j]==1){
00960
00961
00962
00963 ndofe=Tt->give_ndofe (j);
00964
00965 nne=Tt->give_nne (j);
00966 allocv (nne,nodval);
00967 fillv (0.0,nodval);
00968 allocv (ndofe,lv);
00969 fillv (0.0, lv);
00970
00971
00972 sourcenodalvalues (i,nodval);
00973
00974
00975 source_vector (lcid,j,nodval,lv);
00976
00977 cn = new long [ndofe];
00978 Tt->give_code_numbers (j,cn);
00979
00980 locglob (rhs,lv.a,cn,ndofe);
00981
00982 destrv (nodval);
00983 destrv (lv);
00984 delete [] cn;
00985 }
00986
00987 }
00988 }
00989
00990 if (neqs>0){
00991
00992 for (i=0;i<neqs;i++){
00993
00994 j=leqs[i][0];
00995
00996 if (Gtt->leso[j]==1){
00997
00998
00999
01000 ndofe=Tt->give_ndofe (j);
01001
01002 nne=Tt->give_nne (j);
01003 allocv (nne,nodval);
01004 fillv (0.0,nodval);
01005 allocv (ndofe,lv);
01006 fillv (0.0, lv);
01007
01008
01009 k = leqs[i][1];
01010 s = sour[k].giveval (j);
01011 for (k=0;k<nne;k++){
01012 nodval[k]=s;
01013 }
01014
01015
01016 source_vector (lcid,j,nodval,lv);
01017
01018 cn = new long [ndofe];
01019 Tt->give_code_numbers (j,cn);
01020
01021 locglob (rhs,lv.a,cn,ndofe);
01022
01023 destrv (nodval);
01024 destrv (lv);
01025 delete [] cn;
01026 }
01027 }
01028 }
01029 }
01030