00001 #include <math.h>
00002 #include "loadcase.h"
00003 #include "global.h"
00004 #include "globmat.h"
00005 #include "elemswitch.h"
00006 #include "gtopology.h"
00007 #include "matrix.h"
00008 #include "vector.h"
00009 #include "loadn.h"
00010 #include "loadel.h"
00011 #include "element.h"
00012 #include "node.h"
00013 #include "elemhead.h"
00014 #include "intpoints.h"
00015 #include "tablefunct.h"
00016
00017
00018
00019
00020
00021
00022 loadcase::loadcase (void)
00023 {
00024
00025 nln = 0;
00026
00027 nle = 0;
00028
00029 npd = 0;
00030
00031 npt = 0;
00032 tempchang=0;
00033
00034
00035 lon = NULL;
00036
00037 loe = NULL;
00038
00039 pd = NULL;
00040
00041 pt = NULL;
00042
00043 mstress = NULL;
00044
00045 mstrain = NULL;
00046
00047 temp_file = NULL;
00048
00049 pts = NULL;
00050
00051 temp_time = NULL;
00052 }
00053
00054
00055
00056
00057
00058
00059
00060
00061 loadcase::~loadcase (void)
00062 {
00063 long i;
00064
00065 if (lon!=NULL)
00066 delete [] lon;
00067 if (loe!=NULL)
00068 delete [] loe;
00069 if (pd!=NULL)
00070 delete [] pd;
00071 if (pt!=NULL)
00072 delete [] pt;
00073 if (mstress!=NULL)
00074 delete [] mstress;
00075 if (mstrain!=NULL)
00076 delete [] mstrain;
00077 if (temp_file!=NULL)
00078 delete [] temp_file;
00079
00080 if(tempchang == 5){
00081 if (pts!=NULL){
00082 for(i=0;i<ntemp_file;i++)
00083 delete pts[i];
00084 }
00085 }
00086 delete [] pts;
00087
00088 }
00089
00090
00091
00092
00093
00094
00095
00096
00097
00098
00099
00100
00101
00102 void loadcase::read (XFILE *in)
00103 {
00104 long i,j;
00105 XFILE *in_data;
00106
00107
00108 xfscanf (in,"%ld",&nln);
00109 fprintf (stdout,"\n the number of loaded nodes %ld",nln);
00110 lon = new loadn [nln];
00111 for (i=0;i<nln;i++){
00112 lon[i].read (in);
00113 }
00114
00115
00116 xfscanf (in,"%ld",&nle);
00117 fprintf (stdout,"\n the number of loaded elements %ld",nle);
00118 loe = new loadel [nle];
00119 for (i=0;i<nle;i++){
00120 loe[i].read (in);
00121 }
00122
00123
00124 xfscanf (in,"%ld",&npd);
00125 fprintf (stdout,"\n the number of prescribed displacements %ld",npd);
00126 pd = new double [npd];
00127 for (i=0;i<npd;i++){
00128 xfscanf (in,"%lf",pd+i);
00129 }
00130
00131
00132 xfscanf (in,"%ld",&tempchang);
00133 if (tempchang==0)
00134 fprintf (stdout,"\n no temperature changes are described");
00135 switch (tempchang)
00136 {
00137 case 0:
00138 break;
00139 case 1:
00140 case 2:
00141 case 4:
00142 case 5:
00143 if (Mp->temperature < 2)
00144 Mp->temperature=1;
00145 else
00146 {
00147 print_err("type 1 or 2 of temperature load cannot be used together with type 3", __FILE__, __LINE__, __func__);
00148 abort();
00149 }
00150 break;
00151 case 3:
00152 if (Mp->temperature == 1)
00153 {
00154 print_err("type 1 or 2 of temperature load cannot be used together with type 3", __FILE__, __LINE__, __func__);
00155 abort();
00156 }
00157 else
00158 Mp->temperature=2;
00159 break;
00160 default:
00161 print_err("unknown type of temperature load (=%ld) is required", __FILE__, __LINE__, __func__, tempchang);
00162 abort();
00163 }
00164
00165 if (tempchang==1 || tempchang==2){
00166 fprintf (stdout,"\n temperature loading from input file");
00167 npt = Mt->nn;
00168 pt = new double [npt];
00169 for (i=0;i<npt;i++){
00170 xfscanf (in,"%lf",pt+i);
00171 }
00172 }
00173
00174 if (tempchang==3)
00175 fprintf (stdout,"\n temperatures are computing in TRFEL");
00176
00177 if (tempchang==4){
00178 fprintf (stdout,"\n temperature loading - reading temperatures from extra input file");
00179 xfscanf(in, " %a", temp_file);
00180
00181 in_data = xfopen(temp_file,"r");
00182
00183 npt = Mt->nn;
00184 pt = new double [npt];
00185 for (i=0;i<npt;i++){
00186 xfscanf (in_data,"%lf",pt+i);
00187 }
00188 xfclose(in_data);
00189 }
00190
00191 if (tempchang==5){
00192 fprintf (stdout,"\n temperature loading - reading temperatures from extra multiple input files");
00193 xfscanf(in, " %ld",&ntemp_file);
00194 pts = new double* [ntemp_file];
00195 temp_time = new double [ntemp_file];
00196 for (i=0;i<ntemp_file;i++){
00197 xfscanf(in, " %lf",temp_time+i);
00198 xfscanf(in, " %a", temp_file);
00199 in_data = xfopen(temp_file,"r");
00200 npt = Mt->nn;
00201 pts[i] = new double [npt];
00202 for (j=0;j<npt;j++){
00203 xfscanf (in_data,"%lf",&pts[i][j]);
00204 }
00205 xfclose(in_data);
00206 }
00207 }
00208
00209 if (Mp->homog == 3)
00210 {
00211 long ncomp = Mt->give_tncomp(0);
00212 mstress = new double [ncomp];
00213 memset(mstress, 0, sizeof(*mstress)*ncomp);
00214 for (i=0; i<ncomp; i++)
00215 xfscanf(in, "%lf", mstress+i);
00216 }
00217 if (Mp->homog == 4)
00218 {
00219 long ncomp = Mt->give_tncomp(0);
00220 mstrain = new double [ncomp];
00221 memset(mstrain, 0, sizeof(*mstrain)*ncomp);
00222 for (i=0; i<ncomp; i++)
00223 xfscanf(in, "%lf", mstrain+i);
00224 }
00225
00226
00227
00228
00229 }
00230
00231
00232
00233
00234
00235
00236
00237
00238
00239
00240
00241
00242 void loadcase::print (FILE *out)
00243 {
00244 long i;
00245
00246
00247 fprintf(out,"\n## loaded nodes:");
00248 fprintf (out,"\n%ld\n\n",nln);
00249 for (i=0;i<nln;i++){
00250 lon[i].print (out);
00251 }
00252
00253
00254 fprintf(out,"\n## loaded elements:");
00255 fprintf (out,"\n%ld\n\n",nle);
00256 for (i=0;i<nle;i++){
00257 loe[i].print (out,0);
00258 }
00259
00260
00261 fprintf(out,"\n## prescribed displacements:");
00262 fprintf (out,"\n%ld\n\n",npd);
00263 for (i=0;i<npd;i++){
00264 fprintf (out,"\n %lf",pd[i]);
00265 }
00266
00267
00268 fprintf(out,"\n## temperature changes:");
00269 fprintf(out,"\n%ld\n\n",tempchang);
00270
00271 if (tempchang==1 || tempchang==2){
00272 for (i=0;i<npt;i++){
00273 fprintf (out,"\n %lf",pt[i]);
00274 }
00275 }
00276
00277 if (tempchang==4){
00278 fprintf(out, " %s", temp_file);
00279 fprintf(out, "\n");
00280 }
00281
00282 if (tempchang==5){
00283 }
00284
00285 if (Mp->homog == 3)
00286 {
00287 long ncomp = Mt->give_tncomp(0);
00288 for (i=0; i<ncomp; i++)
00289 fprintf(out, "%lf ", mstress[i]);
00290 fprintf(out, "\n");
00291 }
00292 }
00293
00294
00295
00296
00297
00298
00299
00300
00301
00302
00303
00304
00305
00306
00307
00308
00309 void loadcase::assemble (long lcid,double *rhs,double *flv,double scale)
00310 {
00311 long i,j,k,eid,ndofn,ndofe,ndofemn,ne,ncomp,tnip;
00312 ivector cn;
00313 vector r,f,af,anf;
00314 matrix sm;
00315
00316
00317
00318
00319
00320
00321 for (i=0;i<nln;i++){
00322 ndofn=Mt->give_ndofn (lon[i].nid);
00323 if (ndofn<0){
00324 print_err("loaded node %ld is a hanging node",__FILE__,__LINE__,__func__,i++);
00325 }
00326 for (j=0;j<ndofn;j++){
00327 k=Mt->give_dof (lon[i].nid,j);
00328 if (k<=0) continue;
00329 if (k>0) rhs[k-1]+=lon[i].f[j];
00330 }
00331 }
00332
00333
00334
00335
00336
00337 for (i=0;i<nle;i++){
00338
00339 eid=loe[i].eid;
00340 if (Gtm->leso[eid]==1){
00341
00342
00343 ndofe=Mt->give_ndofe (eid);
00344
00345
00346 ndofemn = Gtm->give_ndofe (eid);
00347
00348 reallocv (ndofemn,cn);
00349
00350 Mt->give_code_numbers (eid,cn.a);
00351
00352 if (ndofe != ndofemn){
00353
00354
00355 reallocv (ndofemn,anf);
00356 mtxv (*Mt->elements[i].tmat,loe[i].nf,anf);
00357 locglob (rhs,anf.a,cn.a,ndofemn);
00358 }else{
00359
00360 locglob (rhs,loe[i].nf.a,cn.a,ndofe);
00361 }
00362 }
00363 }
00364
00365 cmulv(scale,rhs,Ndofm);
00366
00367 if (flv)
00368 copyv (rhs,flv,Ndofm);
00369
00370
00371
00372
00373
00374
00375
00376
00377
00378 if ((npd>0) && (Mp->tprob != growing_mech_structure))
00379 {
00380
00381 ne = Mt->ne;
00382 for (i=0;i<ne;i++){
00383 if (Mt->elements[i].prescdispl==1){
00384
00385
00386
00387
00388 ndofe=Mt->give_ndofe (i);
00389 reallocm (ndofe,ndofe,sm);
00390
00391 stiffmat (lcid,i,sm);
00392
00393 reallocv (ndofe,r);
00394 reallocv (ndofe,f);
00395
00396 elprdispl (lcid,i,r.a);
00397
00398 mxv (sm,r,f);
00399 cmulv (-1.0,f);
00400
00401
00402
00403 ndofemn = Gtm->give_ndofe (i);
00404 reallocv (ndofemn,cn);
00405
00406
00407 Mt->give_code_numbers (i,cn.a);
00408
00409 if (ndofe != ndofemn){
00410
00411
00412 reallocv (ndofemn,af);
00413 mtxv (*Mt->elements[i].tmat,f,af);
00414 locglob (rhs,af.a,cn.a,ndofemn);
00415 }else{
00416
00417 locglob (rhs,f.a,cn.a,ndofe);
00418 }
00419
00420 }
00421 }
00422 }
00423
00424
00425 tempercontrib (lcid,rhs,scale);
00426
00427
00428 if (Mp->homog == 3)
00429 {
00430 ncomp = Mt->give_tncomp(0);
00431 j=0;
00432 for (i=Ndofm-ncomp; i<Ndofm; i++)
00433 {
00434 rhs[i] += scale*mstress[j]*Mt->domvol;
00435 if (flv)
00436 flv[i] += scale*mstress[j]*Mt->domvol;
00437 j++;
00438 }
00439 }
00440 if (Mp->homog == 4)
00441 {
00442
00443 ncomp = Mt->give_tncomp(0);
00444
00445 tnip = Mm->tnip;
00446
00447 for (i=0;i<tnip;i++){
00448 for (j=0;j<ncomp;j++){
00449
00450
00451
00452 Mm->ip[i].strain[j]=0.0-mstrain[j];
00453 }
00454
00455 Mm->computenlstresses (i);
00456 }
00457
00458 internal_forces (0,rhs);
00459
00460
00461 }
00462 }
00463
00464
00465
00466
00467
00468
00469
00470
00471
00472
00473
00474
00475
00476 void loadcase::tempercontrib (long lcid,double *rhs,double scale)
00477 {
00478 long i,j;
00479
00480
00481 if (tempchang==1 || tempchang==4){
00482
00483
00484
00485
00486 Mm->est=tempstrain;
00487
00488
00489 intpointval (pt,temperature,1.0);
00490
00491
00492 Mm->temprstrains ();
00493
00494
00495 nodal_eigstrain_forces (lcid,rhs,Mp->time);
00496 }
00497
00498 if (tempchang==2){
00499
00500
00501
00502
00503
00504 Mm->est=tempstrain;
00505
00506
00507 intpointval (pt,temperature,scale);
00508
00509
00510 Mm->temprstrains ();
00511
00512 nodal_eigstrain_forces (lcid,rhs,Mp->time);
00513 }
00514
00515 if (tempchang==3){
00516
00517
00518
00519
00520 Mm->est=tempstrain;
00521
00522
00523 Mm->temprstrains ();
00524
00525
00526 nodal_eigstrain_forces (lcid,rhs,Mp->time);
00527 }
00528
00529
00530 if (tempchang==5){
00531
00532
00533
00534
00535 Mm->est=tempstrain;
00536
00537
00538 npt = Mt->nn;
00539 pt = new double [npt];
00540
00541
00542
00543 tablefunct *tabf1;
00544 tabf1 = new tablefunct;
00545 tabf1->itype = piecewiselin;
00546 tabf1->asize = ntemp_file;
00547 tabf1->x=new double[ntemp_file];
00548 tabf1->y=new double[ntemp_file];
00549 for(i=0;i<Mt->nn;i++){
00550 for(j=0;j<ntemp_file;j++){
00551 tabf1->x[j] = temp_time[j];
00552 tabf1->y[j] = pts[j][i];
00553 }
00554 pt[i] = tabf1->getval(Mp->time);
00555 }
00556
00557
00558
00559 intpointval (pt,temperature,1.0);
00560
00561
00562 Mm->temprstrains ();
00563
00564
00565 nodal_eigstrain_forces (lcid,rhs,Mp->time);
00566 }
00567
00568 }
00569
00570
00571
00572
00573
00574
00575
00576
00577
00578
00579
00580
00581 void loadcase::compute_reactions (long lcid)
00582 {
00583 long i,j,k,ii,nne,ndofn,nn,ne,ndofe,ndofemn,eid;
00584 ivector nod;
00585 vector f,tf;
00586
00587
00588 nn=Mt->nn;
00589
00590 for (i=0;i<nn;i++){
00591 if (Mt->nodes[i].react==1){
00592
00593
00594
00595 ndofn=Mt->give_ndofn (i);
00596 for (k=0;k<ndofn;k++){
00597 Mt->nodes[i].r[k]=0.0;
00598 }
00599 }
00600 }
00601
00602
00603 ne=Mt->ne;
00604
00605 for (i=0;i<ne;i++){
00606 if ((Mt->elements[i].react==1) && (Gtm->leso[i]==1)){
00607
00608
00609
00610
00611 ndofe=Mt->give_ndofe (i);
00612 reallocv (ndofe,f);
00613
00614 elem_internal_forces (i,lcid,f);
00615
00616
00617
00618 ndofemn = Gtm->give_ndofe (i);
00619
00620
00621
00622 nne=Mt->give_nne (i);
00623
00624 if (ndofe != ndofemn){
00625 if (Mp->homog!=3){
00626
00627
00628
00629
00630 vector af;
00631 reallocv (ndofemn,af);
00632
00633 mtxv (*Mt->elements[i].tmat,f,af);
00634 reallocv (ndofemn,f);
00635 copyv (af,f);
00636
00637
00638
00639 nne=Gtm->gelements[i].nmne;
00640 reallocv (nne,nod);
00641 Gtm->gelements[i].give_master_nodes (nod);
00642 }
00643 else{
00644
00645 reallocv (nne,nod);
00646 Mt->give_elemnodes (i,nod);
00647 }
00648 }
00649
00650 if (ndofe == ndofemn){
00651 reallocv (nne,nod);
00652 Mt->give_elemnodes (i,nod);
00653 }
00654
00655
00656
00657 ii=0;
00658
00659
00660 for (j=0;j<nne;j++){
00661
00662
00663 ndofn=Mt->give_ndofn (nod(j));
00664 if (Mt->nodes[nod(j)].react==1){
00665
00666
00667
00668 for (k=0;k<ndofn;k++){
00669 Mt->nodes[nod(j)].r[k]+=f[ii];
00670 ii++;
00671 }
00672 }
00673 else ii+=ndofn;
00674 }
00675
00676 }
00677 }
00678
00679
00680
00681 for (i=0;i<nle;i++){
00682
00683 eid=loe[i].eid;
00684
00685 if ((Mt->elements[eid].react==1) && (Gtm->leso[eid]==1)){
00686
00687
00688
00689 ndofe=Mt->give_ndofe (eid);
00690
00691
00692 ndofemn = Gtm->give_ndofe (eid);
00693
00694 if (ndofe != ndofemn){
00695 if (Mp->homog!=3){
00696
00697
00698
00699
00700 reallocv (ndofemn,f);
00701
00702 mtxv (*Mt->elements[eid].tmat,loe[i].nf,f);
00703
00704
00705
00706 nne=Gtm->gelements[eid].nmne;
00707 reallocv (nne,nod);
00708 Gtm->gelements[eid].give_master_nodes (nod);
00709 }
00710 else{
00711
00712
00713
00714 nne=Mt->give_nne (eid);
00715 reallocv (nne,nod);
00716 Mt->give_elemnodes (eid,nod);
00717 reallocv (ndofe,f);
00718 copyv (loe[i].nf,f);
00719 }
00720 }
00721
00722 if (ndofe == ndofemn){
00723
00724 nne=Mt->give_nne (eid);
00725 reallocv (nne,nod);
00726 Mt->give_elemnodes (eid,nod);
00727 reallocv (ndofe,f);
00728 copyv (loe[i].nf,f);
00729 }
00730
00731
00732
00733 ii=0;
00734
00735
00736 for (j=0;j<nne;j++){
00737
00738
00739 ndofn=Mt->give_ndofn (nod(j));
00740 if (Mt->nodes[nod(j)].react==1){
00741
00742
00743 for (k=0;k<ndofn;k++){
00744 Mt->nodes[nod(j)].r[k]-=f[ii];
00745 ii++;
00746 }
00747 }
00748 else ii+=ndofn;
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