00001 #include "stochdriver.h"
00002 #include "matrix.h"
00003 #include "vector.h"
00004 #include "global.h"
00005 #include "elastisomat.h"
00006 #include "splas1d.h"
00007 #include "j2flow.h"
00008 #include "microM4.h"
00009 #include "chen.h"
00010 #include "crsec2dbar.h"
00011 #include "crsec2dbeam.h"
00012 #include "crsec3dbeam.h"
00013 #include "crsecplstr.h"
00014 #include "crsec3d.h"
00015 #include "loadcase.h"
00016 #include "loadn.h"
00017 #include "intpoints.h"
00018 #include "element.h"
00019 #include "probdesc.h"
00020
00021 stochdriver::stochdriver ()
00022 {
00023
00024 nsmt=0;
00025
00026 nscs=0;
00027
00028 nsnl=0;
00029
00030
00031 nsampl=0;
00032
00033 nstochvar=0;
00034
00035 nprunknowns=0;
00036
00037
00038 mt=NULL; idm=NULL; atm=NULL;
00039
00040 cst=NULL; idcs=NULL; atcs=NULL;
00041
00042 idln=NULL; atln=NULL;
00043
00044
00045 nna=NULL; pnd=NULL;
00046
00047 ena=NULL; ev=NULL;
00048
00049 neigv=0;
00050
00051 avi = NULL; avo = NULL;
00052
00053 datin = NULL;
00054 datout = NULL;
00055
00056
00057 fn = NULL;
00058 }
00059
00060 stochdriver::~stochdriver ()
00061 {
00062 delete [] mt; delete [] idm; delete [] atm;
00063 delete [] cst; delete [] idcs; delete [] atcs;
00064 delete [] idln; delete [] atln;
00065 delete [] nna; delete [] pnd;
00066 delete [] ena; delete [] ev;
00067
00068 delete [] avi; delete [] avo;
00069
00070 delete [] fn;
00071
00072 if (datin!=NULL)
00073 xfclose (datin);
00074 if (datout!=NULL)
00075 fclose (datout);
00076 }
00077
00078
00079
00080
00081
00082
00083
00084
00085 void stochdriver::read (XFILE *in)
00086 {
00087 long i;
00088
00089 if (Mp->stochasticcalc==1 || Mp->stochasticcalc==2){
00090
00091 xfscanf (in," %a",auxfilein);
00092 datin = xfopen (auxfilein,"r");
00093 }
00094 if (Mp->stochasticcalc==3){
00095
00096 fg.read (in);
00097 }
00098
00099
00100 xfscanf (in," %a",auxfileout);
00101 datout = fopen (auxfileout,"w");
00102
00103
00104
00105
00106
00107
00108
00109 xfscanf (in,"%ld",&nsmt);
00110 if (Mespr==1) fprintf (stdout,"\n number of stochastic materials %ld",nsmt);
00111 if (nsmt>0){
00112 mt = new mattype [nsmt];
00113 idm = new long [nsmt];
00114 atm = new atsel[nsmt];
00115 for (i=0;i<nsmt;i++){
00116 xfscanf (in,"%m %ld", &mattype_kwdset, (int*)mt+i, idm+i);
00117 idm[i]--;
00118 atm[i].read(in);
00119 }
00120 }
00121
00122
00123 xfscanf (in,"%ld",&nscs);
00124 if (Mespr==1) fprintf (stdout,"\n number of stochastic cross sections %ld",nscs);
00125 if (nscs>0){
00126 cst = new crsectype [nscs];
00127 idcs = new long [nscs];
00128 atcs = new atsel [nscs];
00129 for (i=0;i<nscs;i++){
00130 xfscanf (in,"%m %ld", &crsectype_kwdset, (int*)cst+i, idcs+i);
00131 idcs[i]--;
00132 atcs[i].read(in);
00133 }
00134 }
00135
00136
00137 xfscanf (in,"%ld",&nsnl);
00138 if (Mespr==1) fprintf (stdout,"\n number of stochastic loaded nodes %ld",nsnl);
00139 if (nsnl>0){
00140 idln = new long [nsnl];
00141 atln = new atsel [nsnl];
00142 for (i=0;i<nsnl;i++){
00143 xfscanf (in,"%ld",idln+i);
00144 idln[i]--;
00145 atln[i].read(in);
00146 }
00147 }
00148
00149
00150
00151
00152
00153
00154 xfscanf (in,"%ld",&npnd);
00155 if (Mespr==1) fprintf (stdout,"\n number of printed nodal displacements %ld",npnd);
00156 if (npnd>0){
00157 nna = new long [npnd];
00158 pnd = new atsel [npnd];
00159 for (i=0;i<npnd;i++){
00160 xfscanf (in,"%ld",nna+i);
00161 nna[i]--;
00162 pnd[i].read (in);
00163 }
00164 }
00165
00166
00167
00168 xfscanf (in,"%ld",&npev);
00169 if (Mespr==1) fprintf (stdout,"\n number of elements with required output %ld",npev);
00170 if (npev>0){
00171 ena = new long [npev];
00172 for (i=0;i<npev;i++){
00173 xfscanf (in,"%ld",ena+i);
00174 ena[i]--;
00175 }
00176 }
00177
00178
00179
00180
00181
00182
00183
00184
00185 if (Mp->stochasticcalc==1 || Mp->stochasticcalc==2){
00186 xfscanf (datin,"%ld %ld",&nsampl,&nstochvar);
00187 }
00188 if (Mp->stochasticcalc==3){
00189 nsampl=fg.tncomb;
00190 nstochvar=fg.nfv;
00191
00192 fprintf (stdout,"\n number of samples %ld",nsampl);
00193 }
00194
00195
00196 avi = new double [nstochvar];
00197
00198
00199 if (Mp->stochasticcalc==1)
00200 readtable (datin);
00201
00202 if (Mp->tprob==eigen_dynamics){
00203 neigv=Mp->eigsol.neigv;
00204 }
00205
00206
00207
00208 compute_nprunknowns ();
00209 fprintf (datout,"%ld\n",nprunknowns);
00210
00211
00212
00213 avo = new double [nprunknowns];
00214
00215
00216 if (Mp->stochasticcalc==1)
00217 allocm (nsampl,nprunknowns,stochtabout);
00218
00219
00220 if (Mp->stochasticcalc==3){
00221 fn = new fuzzynum [nprunknowns];
00222 for (i=0;i<nprunknowns;i++){
00223 fn[i].initiate (fg.nalph,fg.alpha);
00224 }
00225 }
00226
00227 }
00228
00229
00230
00231
00232
00233
00234 void stochdriver::compute_nprunknowns ()
00235 {
00236 long i;
00237
00238 nprunknowns=0;
00239
00240
00241 for (i=0;i<npnd;i++){
00242 nprunknowns+=pnd[i].num;
00243 }
00244 ndispl=nprunknowns;
00245
00246
00247
00248 for (i=0;i<npev;i++){
00249 nprunknowns+=6;
00250 abort ();
00251 }
00252
00253 nelem=nprunknowns-ndispl;
00254
00255 if (Mp->tprob==eigen_dynamics){
00256 nprunknowns*=neigv;
00257 }
00258
00259 }
00260
00261
00262
00263
00264
00265
00266
00267
00268 void stochdriver::readtable (XFILE *in)
00269 {
00270 allocm (nsampl,nstochvar,stochtabin);
00271 readm(in, stochtabin);
00272 }
00273
00274
00275
00276
00277
00278
00279
00280
00281 void stochdriver::writetable ()
00282 {
00283 printm(stochtabout,datout,10,20);
00284 }
00285
00286
00287
00288
00289
00290
00291
00292
00293
00294
00295 void stochdriver::changevalues (long sampleid)
00296 {
00297 assemble_new_values (sampleid);
00298 replace_values ();
00299 }
00300
00301
00302
00303
00304
00305
00306
00307
00308
00309 void stochdriver::assemble_new_values (long sampleid)
00310 {
00311 long i;
00312
00313
00314 if (Mp->stochasticcalc==1){
00315 for (i=0;i<nstochvar;i++){
00316 avi[i]=stochtabin[sampleid][i];
00317 }
00318 }
00319
00320 if (Mp->stochasticcalc==2){
00321 for (i=0;i<nstochvar;i++){
00322 xfscanf (datin,"%le",&avi[i]);
00323 }
00324 }
00325
00326 if (Mp->stochasticcalc==3){
00327 fg.give_new_values (sampleid,avi,Out);
00328 }
00329
00330 }
00331
00332
00333
00334
00335
00336
00337 void stochdriver::replace_values ()
00338 {
00339 long i,j,k,l,m;
00340 vector val;
00341
00342
00343 j=0;
00344 for (i=0;i<nsmt;i++){
00345 k=atm[i].num; m=0;
00346 allocv (k,val);
00347 for (l=j;l<j+k;l++){
00348 val[m]=avi[l];
00349 m++;
00350 }
00351 j+=k;
00352 changematerials (i,val);
00353 destrv (val);
00354 }
00355
00356
00357 for (i=0;i<nscs;i++){
00358 k=atcs[i].num; m=0;
00359 allocv (k,val);
00360 for (l=j;l<j+k;l++){
00361 val[m]=avi[l];
00362 m++;
00363 }
00364 j+=k;
00365 changecrsections (i,val);
00366 destrv (val);
00367 }
00368
00369
00370 for (i=0;i<nsnl;i++){
00371 k=atln[i].num; m=0;
00372 allocv (k,val);
00373 for (l=j;l<j+k;l++){
00374 val[m]=avi[l];
00375 m++;
00376 }
00377 j+=k;
00378 changenodloads (i,val);
00379 destrv (val);
00380 }
00381
00382 }
00383
00384
00385
00386
00387
00388
00389
00390
00391
00392 void stochdriver::changematerials (long id,vector &val)
00393 {
00394 switch (mt[id]){
00395 case elisomat:{
00396 Mm->eliso[idm[id]].changeparam (atm[id],val);
00397 break;
00398 }
00399 case simplas1d:{
00400 Mm->spl1d[idm[id]].changeparam (atm[id],val);
00401 break;
00402 }
00403 case jflow:{
00404 Mm->j2f[idm[id]].changeparam (atm[id],val);
00405 break;
00406 }
00407
00408
00409 case microplaneM4:{
00410 Mm->mpM4[idm[id]].changeparam (atm[id],val);
00411 break;
00412 }
00413
00414 case chenplast:{
00415 Mm->chplast[idm[id]].changeparam (atm[id],val);
00416 break;
00417 }
00418
00419
00420
00421
00422
00423
00424
00425
00426
00427
00428
00429
00430
00431
00432
00433
00434
00435
00436
00437
00438
00439
00440
00441
00442
00443
00444
00445
00446
00447
00448
00449
00450
00451
00452
00453
00454
00455
00456
00457
00458
00459
00460
00461
00462
00463
00464
00465
00466
00467
00468
00469 default:{
00470 fprintf (stderr,"\n\n unknown material type is required in function changematerials (file %s, line %d).\n",__FILE__,__LINE__);
00471 }
00472 }
00473 }
00474
00475 void stochdriver::changecrsections (long id,vector &val)
00476 {
00477 switch (cst[id]){
00478 case nocrosssection:{
00479 break;
00480 }
00481 case csbar2d:{
00482 Mc->cs2dbar[idcs[id]].changeparam (atcs[id],val);
00483 break;
00484 }
00485 case csbeam2d:{
00486 Mc->cs2dbeam[idcs[id]].changeparam (atcs[id],val);
00487 break;
00488 }
00489 case csbeam3d:{
00490 Mc->cs3dbeam[idcs[id]].changeparam (atcs[id],val);
00491 break;
00492 }
00493 case csplanestr:{
00494 Mc->csplstr[idcs[id]].changeparam (atcs[id],val);
00495 break;
00496 }
00497 case cs3dprob:{
00498 Mc->cs3d[idcs[id]].changeparam (atcs[id],val);
00499 break;
00500 }
00501 default:{
00502 fprintf (stderr,"\n\n unknown material type is required in function changecrsections (file %s, line %d).\n",__FILE__,__LINE__);
00503 }
00504 }
00505 }
00506
00507 void stochdriver::changenodloads (long id,vector &val)
00508 {
00509 Mb->lc[0].lon[idln[id]].changeparam (atln[id],val);
00510 }
00511
00512
00513
00514
00515
00516
00517 void stochdriver::extractor ()
00518 {
00519 long i,j,k,ii,ci,nid,nloops;
00520
00521 if (Mp->tprob==eigen_dynamics){
00522 nloops=neigv;
00523 }
00524 else{
00525 nloops=1;
00526 }
00527
00528 ci=0;
00529 for (ii=0;ii<nloops;ii++){
00530
00531 for (i=0;i<npnd;i++){
00532 nid=nna[i];
00533 for (j=0;j<pnd[i].num;j++){
00534 k=Mt->give_dof (nid,pnd[i].atrib[j])-1;
00535 if (k<0){
00536 fprintf (stderr,"\n\n wrong number of DOF in function extractor (file %s, line %d).\n",__FILE__,__LINE__);
00537 }
00538 avo[ci]=Lsrs->lhs[ii*Ndofm+k];
00539 ci++;
00540 }
00541 }
00542
00543
00544
00545
00546
00547
00548
00549
00550
00551
00552
00553
00554
00555
00556
00557
00558
00559
00560
00561
00562
00563
00564
00565
00566
00567
00568 }
00569
00570 }
00571
00572
00573
00574
00575
00576
00577
00578
00579
00580 void stochdriver::save_results (long sampleid)
00581 {
00582 long i;
00583
00584 if (Mp->stochasticcalc==1){
00585 for (i=0;i<nprunknowns;i++){
00586 stochtabout[sampleid][i]=avo[i];
00587 }
00588 }
00589 if (Mp->stochasticcalc==2){
00590 for (i=0;i<nprunknowns;i++){
00591 fprintf (datout,"%le ",avo[i]);
00592 }
00593 fprintf (datout,"\n");
00594 }
00595 if (Mp->stochasticcalc==3){
00596 fg.save_values (fn,nprunknowns,avo);
00597 }
00598
00599 }
00600
00601
00602
00603
00604
00605
00606
00607 void stochdriver::diagpostproc ()
00608 {
00609 long i,j,nc,nr,minnr;
00610 double x,y,minx,dx;
00611 char *oname,fname[1000];
00612 char *path,*name,*suffix;
00613 FILE *in,*out;
00614 gfunct *gf;
00615
00616
00617 nc = Outdm->odiag[0].npun;
00618
00619 if (nc!=2){
00620 fprintf (stderr,"\n\n unsupported number of columns in function diagpostproc (file %s, line %d)\n",__FILE__,__LINE__);
00621 abort ();
00622 }
00623
00624
00625 oname = new char [1000];
00626 for (i=0;i<1000;i++){
00627 oname[i] = Outdm->outdiagfn[i];
00628 }
00629
00630 gf = new gfunct [nsampl];
00631
00632
00633
00634 minx=1.0e8;
00635 minnr=10000000;
00636 for (i=0;i<nsampl;i++){
00637
00638 filename_decomposition (oname,path,name,suffix);
00639
00640 sprintf(fname, "%s.%ld.dat",name,i+1);
00641 in=fopen (fname,"r");
00642
00643 delete [] path;
00644 delete [] name;
00645 delete [] suffix;
00646
00647
00648 nr=0;
00649 do{
00650 fscanf (in,"%le %le",&x,&x);
00651 nr++;
00652
00653
00654
00655 j=feof(in);
00656 }while (j==0);
00657
00658 nr--;
00659 if (minnr>nr) minnr=nr;
00660
00661 rewind (in);
00662 gf[i].init_tab (nr);
00663
00664 for (j=0;j<nr;j++){
00665 fscanf (in,"%le %le",&x,&y);
00666 gf[i].tabf->x[j]=fabs(x);
00667 gf[i].tabf->y[j]=fabs(y);
00668 }
00669
00670 if (fabs(x)<minx) minx=fabs(x);
00671
00672 fclose (in);
00673 }
00674
00675 printf ("\n minnr %ld minx %le\n",minnr,minx);
00676
00677
00678
00679
00680 out = fopen ("fuzzy.dat","w");
00681 dx=minx/minnr;
00682 x=1.0e-8;
00683 for (i=0;i<minnr;i++){
00684 fprintf (out,"%le ",x);
00685 for (j=0;j<nsampl;j++){
00686 fprintf (out,"%le ",gf[j].getval(x));
00687 }
00688 fprintf (out,"\n");
00689 x+=dx;
00690 }
00691 fclose (out);
00692
00693
00694
00695
00696
00697
00698
00699
00700
00701 }
00702
00703 void stochdriver::update_auxparam ()
00704 {
00705 if (Mp->stochasticcalc==3){
00706
00707
00708 fg.actcomb++;
00709 if (fg.actcomb==fg.ncomb){
00710 fg.actalph++;
00711 fg.actcomb=0;
00712 }
00713
00714 }
00715
00716 }
00717
00718
00719
00720
00721
00722
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
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