00001 #include "fdsolver.h"
00002 #include "edsolver.h"
00003 #include "global.h"
00004 #include "globmat.h"
00005 #include "gmatrix.h"
00006 #include "loadcase.h"
00007 #include "dloadcase.h"
00008 #include "mechprint.h"
00009 #include "elemswitch.h"
00010
00011
00012
00013 #include "intpoints.h"
00014 #include "element.h"
00015
00016
00017 void solve_forced_dynamics ()
00018 {
00019 long i;
00020
00021 for (i=0;i<Lsrs->nlc;i++){
00022 switch (Mp->tforvib){
00023 case newmark:{
00024 newmark_method (i);
00025 break;
00026 }
00027 case findiff:{
00028 difference_method (i);
00029 break;
00030 }
00031 case modal_analysis:{
00032 response_spectrum_method (i);
00033 break;
00034 }
00035 default:{
00036 fprintf (stderr,"\n\n unknown solver of forced vibration is required");
00037 fprintf (stderr,"\n in function solve_forced_dynamics () (file %s, line %d).\n",__FILE__,__LINE__);
00038 }
00039 }
00040 }
00041
00042 }
00043
00044
00045
00046
00047
00048
00049
00050
00051
00052 void newmark_method (long lcid)
00053 {
00054 long i,j,n;
00055 double alpha,delta,zero,time,dt,end_time;
00056 double *a,*v,*p,*dd,*vv,*lhs,*rhs,*fs;
00057
00058
00059 n = Ndofm;
00060
00061 zero=Mp->zero;
00062
00063 alpha=Mp->alphafvn;
00064
00065 delta=Mp->deltafvn;
00066
00067
00068 time=Mp->timecon.starttime ();
00069
00070 end_time = Mp->timecon.endtime ();
00071
00072 dt = Mp->timecon.initialtimeincr ();
00073
00074
00075
00076 lhs = Lsrs->give_lhs (lcid);
00077
00078
00079 v = Lsrs->give_tdlhs (lcid);
00080
00081 rhs = Lsrs->give_rhs (2*lcid);
00082
00083 fs = Lsrs->give_rhs (2*lcid+1);
00084
00085
00086 a = Lsrs->give_stdlhs (lcid);
00087 nullv (a,n);
00088
00089 dd = new double [n];
00090
00091 vv = new double [n];
00092
00093 p = new double [n];
00094
00095
00096 nullv (fs+lcid*n, n);
00097 Mb->lc[lcid].assemble (lcid,fs,NULL,1.0);
00098
00099
00100 i=0;
00101 print_init(-1, "wt");
00102 print_step(lcid, i, Mp->time, rhs);
00103 print_flush();
00104
00105
00106
00107 Mp->time = Mp->timecon.newtime ();
00108
00109 dt = Mp->timecon.actualbacktimeincr ();
00110
00111
00112
00113 do{
00114
00115
00116 stiffness_matrix (lcid);
00117
00118 mass_matrix (lcid);
00119
00120 damping_matrix ();
00121
00122 if (Mespr==1) fprintf (stdout,"\n time %f",Mp->time);
00123
00124
00125 nullv (rhs+lcid*n, n);
00126 Mb->dlc[lcid].assemble (lcid,rhs,NULL,Ndofm,Mp->time);
00127
00128
00129 addv (rhs,fs,n);
00130
00131
00132 cmulv (dt*dt*alpha,rhs,n);
00133
00134
00135 for (j=0;j<n;j++){
00136 dd[j] = lhs[j] + dt*v[j] + dt*dt*(0.5-alpha)*a[j];
00137 vv[j] = v[j] + dt*(1.0-delta)*a[j];
00138 }
00139
00140
00141 Mmat->gmxv (dd,p);
00142
00143
00144 addv (rhs,p,n);
00145
00146
00147 Dmat->gmxv (dd,p);
00148 cmulv (dt*delta,p,n);
00149 addv (rhs,p,n);
00150
00151
00152 Dmat->gmxv (vv,p);
00153 cmulv (-1.0*dt*dt*alpha,p,n);
00154 addv (rhs,p,n);
00155
00156
00157
00158
00159 Smat->scalgm (dt*dt*alpha);
00160 Smat->addgm (1.0,*Mmat);
00161 Smat->addgm (dt*delta,*Dmat);
00162
00163
00164
00165
00166 Mp->ssle->solve_system (Gtm,Smat,lhs,rhs,Out);
00167
00168 for (j=0;j<n;j++){
00169 v[j] = vv[j] + delta/dt/alpha*(lhs[j]-dd[j]);
00170 a[j] = 1.0/dt/dt/alpha*(lhs[j]-dd[j]);
00171 }
00172
00173 compute_req_val (lcid);
00174 print_step(lcid, i, Mp->time, rhs);
00175 print_flush();
00176
00177
00178 Mp->time = Mp->timecon.newtime ();
00179 dt = Mp->timecon.actualbacktimeincr ();
00180 i++;
00181
00182 }while(Mp->time<end_time);
00183
00184 print_close();
00185 delete [] p;
00186 delete [] dd;
00187 delete [] vv;
00188 }
00189
00190
00191
00192
00193
00194
00195
00196
00197
00198
00199 void difference_method (long lcid)
00200 {
00201 long i,j,n;
00202 double zero,time,dt,end_time;
00203 double *a,*v,*p,*da,*dp,*lhs,*rhs,*fs;
00204
00205
00206 n = Ndofm;
00207
00208 zero=Mp->zero;
00209
00210
00211 time=Mp->timecon.starttime ();
00212
00213 end_time = Mp->timecon.endtime ();
00214
00215 dt = Mp->timecon.initialtimeincr ();
00216
00217
00218
00219 lhs = Lsrs->give_lhs (lcid);
00220
00221
00222 v = Lsrs->give_tdlhs (lcid);
00223
00224 a = Lsrs->give_stdlhs (lcid);
00225
00226 rhs = Lsrs->give_rhs (2*lcid);
00227
00228 fs = Lsrs->give_rhs (2*lcid+1);
00229
00230
00231
00232 dp = new double [n];
00233
00234 da = new double [n];
00235
00236 p = new double [n];
00237
00238
00239 nullv (fs+lcid*n, n);
00240 Mb->lc[lcid].assemble (lcid,fs,NULL,1.0);
00241
00242
00243 i=0;
00244 print_init(-1, "wt");
00245 print_step(lcid, i, Mp->time, rhs);
00246 print_flush();
00247
00248
00249
00250 for (j=0;j<n;j++){
00251 da[j] = lhs[j];
00252 dp[j] = da[j] - dt*v[j];
00253 }
00254
00255
00256
00257 do{
00258
00259
00260 stiffness_matrix (lcid);
00261
00262 mass_matrix (lcid);
00263
00264 damping_matrix ();
00265
00266 if (Mespr==1) fprintf (stdout,"\n time %f",Mp->time);
00267
00268
00269 nullv (rhs+lcid*n, n);
00270 Mb->dlc[lcid].assemble (lcid,rhs,NULL,Ndofm,Mp->time);
00271
00272
00273 addv (rhs,fs,n);
00274
00275
00276 Smat->gmxv (da,p);
00277 cmulv (-1.0,p,n);
00278 addv (rhs,p,n);
00279
00280
00281 Mmat->gmxv (da,p);
00282 cmulv (2.0/dt/dt,p,n);
00283 addv (rhs,p,n);
00284
00285
00286 Mmat->gmxv (dp,p);
00287 cmulv (-1.0/dt/dt,p,n);
00288 addv (rhs,p,n);
00289
00290
00291 Dmat->gmxv (dp,p);
00292 cmulv (1.0/2.0/dt,p,n);
00293 addv (rhs,p,n);
00294
00295
00296
00297
00298 Mmat->scalgm (1.0/dt/dt);
00299 Mmat->addgm (1.0/2.0/dt,*Dmat);
00300
00301
00302
00303 Mp->ssle->solve_system (Gtm,Mmat,lhs,rhs,Out);
00304
00305 for (j=0;j<n;j++){
00306 v[j]=(lhs[j]-dp[j])/2.0/dt;
00307 a[j]=(lhs[j]-2.0*da[j]+dp[j])/dt/dt;
00308
00309 dp[j]=da[j];
00310 da[j]=lhs[j];
00311 }
00312
00313
00314 Mp->time = Mp->timecon.newtime ();
00315 dt = Mp->timecon.actualbacktimeincr ();
00316 i++;
00317
00318 compute_req_val (lcid);
00319 print_step(lcid, i, Mp->time, rhs);
00320 print_flush();
00321
00322
00323 }while(Mp->time<end_time);
00324
00325 print_close();
00326 delete [] da;
00327 delete [] dp;
00328 delete [] p;
00329 }
00330
00331
00332
00333
00334
00335
00336
00337
00338
00339
00340
00341
00342
00343 void verlet_method (long lcid)
00344 {
00345 long i,n,ni;
00346 double mm,zero,dt,end_time;
00347 double *m,*dp,*da,*dn,*fp,*fi;
00348
00349
00350 n = Ndofm;
00351
00352 ni = 0;
00353
00354 zero = Mp->zero;
00355
00356
00357 dp = new double [n];
00358
00359 dn = new double [n];
00360
00361 da = Lsrs->give_lhs (lcid);
00362
00363 fp = Lsrs->give_rhs (lcid);
00364
00365 fi = new double [n];
00366
00367 m = new double [n];
00368
00369
00370 Mp->time=Mp->timecon.starttime ();
00371
00372 end_time = Mp->timecon.endtime ();
00373
00374 dt = Mp->timecon.initialtimeincr ();
00375
00376
00377 for (i=0;i<n;i++){
00378 mm=Mmat->give_entry (i,i);
00379 if (mm<zero){
00380
00381 }
00382 else{
00383 m[i]=1.0/mm;
00384 }
00385 }
00386
00387
00388
00389
00390 mefel_right_hand_side (i,fp);
00391
00392
00393 internal_forces (lcid,fi);
00394
00395 for (i=0;i<n;i++){
00396
00397 da[i]=Lsrs->lhsi[i];
00398
00399 dn[i]=Lsrs->tdlhsi[i];
00400
00401 dp[i]=da[i]-dt*dn[i]+dt*dt/2.0*m[i]*(fp[i]+fi[i]);
00402 }
00403
00404
00405 print_init(-1, "wt");
00406
00407 print_flush();
00408
00409
00410 do{
00411 fprintf (stdout,"\n time %e",Mp->time);
00412
00413
00414 mefel_right_hand_side (i,fp);
00415
00416
00417 internal_forces (lcid,fi);
00418
00419 for (i=0;i<n;i++){
00420 dn[i]=2.0*da[i]-dp[i]+dt*dt*m[i]*(fi[i]+fp[i]);
00421 dp[i]=da[i];
00422 da[i]=dn[i];
00423 }
00424
00425
00426
00427 print_flush();
00428
00429
00430 Mp->time = Mp->timecon.newtime ();
00431 ni++;
00432
00433
00434
00435 }while(Mp->time<end_time);
00436
00437 print_close ();
00438
00439 }
00440
00441
00442
00443
00444
00445
00446
00447
00448 void response_spectrum_method (long lcid)
00449 {
00450 long i,j,n,m;
00451 double coeff,per;
00452 double *v,*auxv;
00453
00454
00455 n=Ndofm;
00456
00457 m=Mp->eigsol.neigv;
00458
00459
00460 solve_eigen_dynamics (Lsrs->lhs+n,Lsrs->w);
00461
00462 v = new double [n];
00463 auxv = new double [n];
00464
00465 nullv (Lsrs->lhs,n);
00466
00467 for (i=0;i<m;i++){
00468
00469 per=2.0*3.14159265358979/sqrt(Lsrs->w[i]);
00470
00471 nullv (v,n);
00472
00473
00474 Mb->dlc[0].stool.assemble (v,per);
00475
00476
00477 Mmat->gmxv (v,auxv);
00478
00479
00480 coeff = ss (auxv,Lsrs->lhs+(i+1)*n,n)/(Lsrs->w[i]);
00481
00482
00483
00484 for (j=0;j<n;j++){
00485 Lsrs->lhs[j]-=coeff*Lsrs->lhs[(i+1)*n+j];
00486 }
00487
00488 }
00489
00490
00491 for (j=0;j<n;j++){
00492 Lsrs->lhs[j]=Lsrs->lhs[n+j];
00493 }
00494
00495
00496
00497
00498
00499
00500
00501
00502 compute_ipstresses (0);
00503
00504
00505
00506
00507
00508
00509
00510
00511
00512
00513
00514
00515 delete [] auxv;
00516 delete [] v;
00517 }