00001 #include <math.h>
00002 #include "dloadcase.h"
00003 #include "loadcase.h"
00004 #include "global.h"
00005 #include "gtopology.h"
00006 #include "matrix.h"
00007 #include "vector.h"
00008 #include "globmat.h"
00009 #include "dloadn.h"
00010 #include "dloadpd.h"
00011 #include "element.h"
00012 #include "node.h"
00013 #include "elemswitch.h"
00014
00015
00016
00017
00018
00019
00020
00021
00022 dloadcase::dloadcase (void)
00023 {
00024
00025 nslc=0;
00026
00027 nln = 0;
00028
00029 nle = 0;
00030
00031 npd = 0;
00032
00033 lon = NULL; pd = NULL;
00034
00035 slc = NULL;
00036
00037 gf=NULL;
00038 }
00039
00040
00041
00042
00043
00044
00045
00046
00047 dloadcase::~dloadcase (void)
00048 {
00049 delete [] lon;
00050 delete [] pd;
00051 delete [] slc;
00052 delete [] gf;
00053 }
00054
00055
00056
00057
00058
00059
00060
00061
00062
00063
00064 void dloadcase::read (XFILE *in)
00065 {
00066 long i;
00067
00068
00069 xfscanf (in,"%k%m","type_of_time_load",&dynload_kwdset,(int*)&tdl);
00070
00071 switch (tdl){
00072 case timeindload:{
00073
00074
00075
00076 xfscanf (in,"%ld",&nslc);
00077
00078 slc = new loadcase [nslc];
00079 gf = new gfunct [nslc];
00080
00081 for (i=0;i<nslc;i++){
00082
00083 slc[i].read (in);
00084
00085 gf[i].read (in);
00086 }
00087 break;
00088 }
00089 case seismicload:{
00090
00091 stool.read (in);
00092 break;
00093 }
00094 case responsespectrum:{
00095
00096 stool.read (in);
00097 break;
00098 }
00099 case timedepload:{
00100
00101
00102
00103
00104
00105 xfscanf (in,"%ld",&nln);
00106 lon = new dloadn [nln];
00107 for (i=0;i<nln;i++){
00108 lon[i].read (in);
00109 }
00110
00111
00112
00113
00114
00115
00116
00117
00118
00119
00120 xfscanf (in,"%ld",&npd);
00121 pd = new dloadpd [npd];
00122 for (i=0;i<npd;i++){
00123 pd[i].read (in);
00124 }
00125 break;
00126 }
00127
00128 default:{
00129 print_err("unknown type of load case is required", __FILE__, __LINE__, __func__);
00130 }
00131 }
00132
00133 }
00134
00135
00136
00137
00138
00139
00140
00141
00142
00143
00144 void dloadcase::print (FILE *out)
00145 {
00146 long i;
00147
00148
00149 fprintf (out,"\n ##type of time-dependent load:");
00150 fprintf (out,"\n %d\n",(int)tdl);
00151
00152 switch (tdl){
00153 case timeindload:{
00154
00155
00156 fprintf (out,"\n #time-independent load case:");
00157 fprintf (out,"\n #number of subloadcases:");
00158 fprintf (out,"\n\n %ld",nslc);
00159
00160 for (i=0;i<nslc;i++){
00161
00162 slc[i].print (out);
00163
00164 gf[i].print (out);
00165 }
00166 break;
00167 }
00168 case seismicload:{
00169
00170 fprintf (out,"\n #seismic load:");
00171
00172 break;
00173 }
00174 case responsespectrum:{
00175
00176 fprintf (out,"\n #response spectrum:");
00177
00178 break;
00179 }
00180 case timedepload:{
00181
00182
00183
00184
00185
00186 fprintf (out,"\n #time-dependent load case:");
00187 fprintf (out,"\n #loaded nodes:");
00188 fprintf (out,"\n %ld",nln);
00189 for (i=0;i<nln;i++){
00190 lon[i].print (out);
00191 }
00192
00193
00194
00195
00196
00197
00198
00199
00200
00201
00202 fprintf (out,"\n #prescribed displacements:");
00203 fprintf (out,"\n\n %ld",npd);
00204 for (i=0;i<npd;i++){
00205 pd[i].print (out);
00206 }
00207 break;
00208 }
00209
00210 default:{
00211 print_err("unknown type of load case is required", __FILE__, __LINE__, __func__);
00212 }
00213 }
00214
00215 }
00216
00217
00218
00219
00220
00221
00222
00223
00224
00225
00226
00227
00228
00229
00230
00231
00232
00233 void dloadcase::assemble (long lcid,double *rhs,double *flv,long n,double t)
00234 {
00235 long i,j,k,ndofn,ndofe,ndofemn,slcpd;
00236 ivector cn;
00237 vector r,f,af;
00238 matrix sm;
00239 double a,*aux,*auxf;
00240
00241
00242
00243 switch (tdl){
00244 case timeindload:{
00245
00246
00247
00248 aux = new double [n];
00249 auxf = new double [n];
00250 slcpd = 0;
00251 for (i=0;i<nslc;i++){
00252 nullv (aux,n);
00253 nullv (auxf,n);
00254
00255 if (slc[i].npd)
00256 slcpd++;
00257
00258 a=gf[i].getval(t);
00259
00260
00261 slc[i].assemble (0,aux,auxf,a);
00262
00263 for (j=0;j<n;j++){
00264 rhs[j]+=aux[j];
00265 }
00266
00267 if (flv)
00268 {
00269 for (j=0;j<n;j++)
00270 flv[j]+=auxf[j];
00271 }
00272 }
00273
00274
00275
00276
00277
00278 if (slcpd>0){
00279
00280 long ne = Mt->ne;
00281
00282 for (i=0;i<ne;i++){
00283 if (Mt->elements[i].prescdispl==1){
00284
00285
00286
00287
00288 ndofe=Mt->give_ndofe (i);
00289 reallocm (ndofe,ndofe,sm);
00290
00291 stiffmat (lcid,i,sm);
00292
00293 reallocv (ndofe,r);
00294 reallocv (ndofe,f);
00295
00296
00297
00298 elprdispl (lcid,i,r.a);
00299
00300 mxv (sm,r,f);
00301 cmulv (-1.0,f);
00302
00303
00304
00305 ndofemn = Gtm->give_ndofe (i);
00306 reallocv (ndofemn,cn);
00307
00308
00309 Mt->give_code_numbers (i,cn.a);
00310
00311 if (ndofe != ndofemn){
00312
00313
00314 reallocv (ndofemn,af);
00315 mtxv (*Mt->elements[i].tmat,f,af);
00316 locglob (rhs,af.a,cn.a,ndofemn);
00317 }else{
00318
00319 locglob (rhs,f.a,cn.a,ndofe);
00320 }
00321
00322 }
00323 }
00324 }
00325
00326 delete [] aux;
00327 delete [] auxf;
00328
00329 break;
00330 }
00331 case seismicload:{
00332
00333
00334
00335
00336 aux = new double [n];
00337
00338 stool.assemble (rhs,t);
00339
00340 Mmat->gmxv (rhs,aux);
00341
00342 copyv (aux,rhs,n);
00343
00344
00345 if (flv)
00346 copyv (aux,flv,n);
00347
00348 delete [] aux;
00349
00350 break;
00351 }
00352 case responsespectrum:{
00353
00354 break;
00355 }
00356 case timedepload:{
00357
00358
00359
00360
00361
00362 for (i=0;i<nln;i++){
00363
00364 ndofn=Mt->give_ndofn (lon[i].idn);
00365 if (ndofn<0){
00366 print_err("loaded node %ld is a hanging node",__FILE__,__LINE__,__func__,i++);
00367 }
00368
00369 for (j=0;j<ndofn;j++){
00370 k=Mt->give_dof (lon[i].idn,j);
00371 if (k<=0) continue;
00372 if (k>0) rhs[k-1]+=lon[i].getval(t,j);
00373 }
00374 }
00375
00376
00377
00378
00379
00380
00381
00382
00383
00384
00385
00386
00387
00388
00389
00390
00391
00392
00393 if (flv)
00394 copyv(rhs,flv,n);
00395
00396
00397
00398
00399 if (npd>0){
00400
00401 long ne = Mt->ne;
00402
00403 for (i=0;i<ne;i++){
00404 if (Mt->elements[i].prescdispl==1){
00405
00406
00407
00408
00409 ndofe=Mt->give_ndofe (i);
00410 reallocm (ndofe,ndofe,sm);
00411
00412 stiffmat (lcid,i,sm);
00413
00414 reallocv (ndofe,r);
00415 reallocv (ndofe,f);
00416
00417
00418 elprdispl (lcid,i,r.a);
00419
00420 mxv (sm,r,f);
00421 cmulv (-1.0,f);
00422
00423
00424
00425 ndofemn = Gtm->give_ndofe (i);
00426 reallocv (ndofemn,cn);
00427
00428
00429 Mt->give_code_numbers (i,cn.a);
00430
00431 if (ndofe != ndofemn){
00432
00433
00434 reallocv (ndofemn,af);
00435 mtxv (*Mt->elements[i].tmat,f,af);
00436 locglob (rhs,af.a,cn.a,ndofemn);
00437 }else{
00438
00439 locglob (rhs,f.a,cn.a,ndofe);
00440 }
00441
00442 }
00443 }
00444 }
00445
00446 break;
00447 }
00448
00449 default:{
00450 print_err("unknown type of load case is required",__FILE__,__LINE__,__func__);
00451 }
00452 }
00453
00454 }
00455
00456
00457
00458
00459
00460
00461
00462
00463
00464
00465
00466
00467
00468 void dloadcase::assemble (double *rhs, double *lhs, double *flv)
00469 {
00470 long i,j,k,ndofn;
00471
00472
00473
00474
00475
00476
00477 for (i=0;i<nln;i++){
00478 ndofn=Mt->give_ndofn (lon[i].idn);
00479 for (j=0;j<ndofn;j++){
00480 k=Mt->give_dof (lon[i].idn,j);
00481 if (k<=0) continue;
00482 if (k>0) rhs[k-1]+=lon[i].getval(lhs[k-1], j);
00483 }
00484 }
00485
00486 if (flv)
00487 copyv(rhs,flv,Ndofm);
00488 }
00489
00490
00491
00492
00493
00494
00495
00496
00497
00498
00499
00500
00501 void dloadcase::compute_reactions (long lcid)
00502 {
00503 long i,j,k,ii,nne,ndofn,ne,ndofe,ndofemn;
00504 vector f,tf,af;
00505 ivector nod;
00506
00507 for (i=0;i<Mt->nn;i++){
00508 if (Mt->nodes[i].react==1){
00509 ndofn=Mt->give_ndofn (i);
00510 for (k=0;k<ndofn;k++){
00511 Mt->nodes[i].r[k]=0.0;
00512 }
00513 }
00514 }
00515
00516
00517 ne=Mt->ne;
00518
00519 for (i=0;i<ne;i++){
00520 if ((Mt->elements[i].react==1) && (Gtm->leso[i]==1)){
00521
00522
00523
00524
00525 ndofe=Mt->give_ndofe (i);
00526 reallocv(ndofe,f);
00527
00528 elem_internal_forces (i,lcid,f);
00529
00530
00531
00532 ndofemn = Gtm->give_ndofe (i);
00533
00534
00535
00536 nne=Mt->give_nne (i);
00537
00538 if (ndofe != ndofemn){
00539
00540
00541
00542 reallocv (ndofemn,af);
00543
00544 mtxv (*Mt->elements[i].tmat,f,af);
00545 reallocv (ndofemn,f);
00546 copyv (af,f);
00547
00548
00549
00550 nne=Gtm->gelements[i].nmne;
00551 reallocv (nne,nod);
00552 Gtm->gelements[i].give_master_nodes (nod);
00553
00554 }
00555
00556 if (ndofe == ndofemn){
00557 reallocv (nne,nod);
00558 Mt->give_elemnodes (i,nod);
00559 }
00560
00561
00562
00563
00564 ii=0;
00565
00566
00567 for (j=0;j<nne;j++){
00568
00569
00570 ndofn=Mt->give_ndofn (nod(j));
00571 if (Mt->nodes[nod(j)].react==1){
00572
00573
00574
00575 for (k=0;k<ndofn;k++){
00576 Mt->nodes[nod(j)].r[k]+=f[ii];
00577 ii++;
00578 }
00579 }
00580 else ii+=ndofn;
00581 }
00582
00583 }
00584 }
00585
00586 }
00587
00588
00589
00590
00591
00592
00593
00594
00595
00596
00597
00598
00599
00600 double dloadcase::get_pd(double time, long dof)
00601 {
00602 double ret = 0.0;
00603 double a;
00604 long i;
00605
00606 switch (tdl){
00607 case timeindload:{
00608 for (i=0; i<nslc; i++)
00609 {
00610
00611 a=gf[i].getval(time);
00612
00613 ret += slc[i].pd[0-dof-1]*a;
00614 }
00615 break;
00616 }
00617 case timedepload:{
00618 ret = pd[0-dof-1].getval(Mp->time);
00619 break;
00620 }
00621 default:{
00622 print_err("this type of load case is not supported", __FILE__, __LINE__, __func__);
00623 abort();
00624 }
00625 }
00626 return ret;
00627 }
00628
00629
00630
00631
00632
00633
00634
00635
00636
00637
00638
00639
00640
00641
00642 void dloadcase::tempercontrib (long lcid,double *rhs,long n,double t)
00643 {
00644 long i;
00645 double a;
00646 double *aux;
00647
00648 switch (tdl){
00649 case timeindload:{
00650 aux = new double [n];
00651 for (i=0;i<nslc;i++){
00652 nullv (aux,n);
00653
00654 a=gf[i].getval(t);
00655
00656
00657 slc[i].tempercontrib (lcid,aux,a);
00658 addv (rhs,aux,n);
00659 }
00660 delete [] aux;
00661 break;
00662 }
00663 default:
00664 print_err("unknown type of load case is required", __FILE__, __LINE__, __func__);
00665 }
00666
00667 }