00001 #include "ndsolver.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 void solve_nonlinear_dynamics ()
00012 {
00013 switch (Mp->tforvib){
00014 case newmark:{
00015 nonlin_newmark_method ();
00016 break;
00017 }
00018 default:{
00019 fprintf (stderr,"\n\n unknown solver of nonlinear dynamics is required");
00020 fprintf (stderr,"\n in function solve_nonlinear_dynamics () (file %s, line %d).\n",__FILE__,__LINE__);
00021 }
00022 }
00023 }
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033
00034 void nonlin_newmark_method ()
00035 {
00036 long i,j,k,n,ini,lcid;
00037 double beta,gamma,zero,time,dt,dtc,end_time,err,norres;
00038 double *a,*v,*p,*dd,*vv,*lhs,*rhs,*fs,*brhs;
00039
00040
00041 n = Ndofm;
00042
00043 beta=Mp->alphafvn;
00044
00045 gamma=Mp->deltafvn;
00046
00047
00048
00049 ini = Mp->nlman->niilnr;
00050
00051 zero = Mp->zero;
00052
00053 err = Mp->nlman->errnr;
00054
00055
00056
00057 time=Mp->timecon.starttime ();
00058
00059 end_time = Mp->timecon.endtime ();
00060
00061 dt = Mp->timecon.initialtimeincr ();
00062
00063 dtc = 0.0;
00064
00065
00066
00067 lcid=0;
00068
00069
00070
00071 lhs = Lsrs->give_lhs (lcid);
00072
00073
00074 v = Lsrs->give_tdlhs (lcid);
00075
00076 rhs = Lsrs->give_rhs (2*lcid);
00077
00078 fs = Lsrs->give_rhs (2*lcid+1);
00079
00080 brhs = new double [n];
00081
00082
00083 a = Lsrs->give_stdlhs (lcid);
00084 nullv (a,n);
00085
00086 dd = new double [n];
00087
00088 vv = new double [n];
00089
00090 p = new double [n];
00091
00092
00093 nullv (fs+lcid*n, n);
00094 Mb->lc[lcid].assemble (lcid,fs,NULL,1.0);
00095
00096
00097 i=0;
00098 print_init(-1, "wt");
00099 print_step(lcid, i, Mp->time, rhs);
00100 print_flush();
00101
00102
00103
00104 Mp->time = Mp->timecon.newtime ();
00105
00106 dt = Mp->timecon.actualbacktimeincr ();
00107
00108
00109
00110 do{
00111
00112
00113 stiffness_matrix (lcid);
00114
00115 mass_matrix (lcid);
00116
00117 damping_matrix ();
00118
00119 if (Mespr==1) fprintf (stdout,"\n time %f",Mp->time);
00120
00121
00122 nullv (rhs+lcid*n, n);
00123 Mb->dlc[lcid].assemble (lcid,rhs,NULL,n,Mp->time);
00124
00125
00126 addv (rhs,fs,n);
00127
00128
00129 cmulv (dt*dt*beta,rhs,n);
00130
00131 copyv (rhs,brhs,n);
00132
00133
00134
00135 for (j=0;j<n;j++){
00136 dd[j] = lhs[j] + dt*v[j] + dt*dt*(0.5-beta)*a[j];
00137 vv[j] = v[j] + dt*(1.0-gamma)*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*gamma,p,n);
00149 addv (rhs,p,n);
00150
00151
00152 Dmat->gmxv (vv,p);
00153 cmulv (-1.0*dt*dt*beta,p,n);
00154 addv (rhs,p,n);
00155
00156
00157
00158
00159 Smat->scalgm (dt*dt*beta);
00160 Smat->addgm (1.0,*Mmat);
00161 Smat->addgm (dt*gamma,*Dmat);
00162
00163
00164
00165 Mp->ssle->solve_system (Gtm,Smat,lhs,rhs,Out);
00166
00167 for (j=0;j<n;j++){
00168 a[j] = 1.0/dt/dt/beta*(lhs[j]-dd[j]);
00169 v[j] = vv[j] + gamma*dt*a[j];
00170 }
00171
00172
00173
00174 copyv (brhs,rhs,n);
00175
00176
00177
00178
00179 Mmat->gmxv (a,p);
00180 cmulv (-1.0*dt*dt*beta,p,n);
00181 addv (rhs,p,n);
00182
00183
00184
00185 Dmat->gmxv (v,p);
00186 cmulv (-1.0*dt*dt*beta,p,n);
00187 addv (rhs,p,n);
00188
00189
00190 internal_forces (lcid,p);
00191 cmulv (-1.0*dt*dt*beta,p,n);
00192 addv (rhs,p,n);
00193
00194
00195 norres = ss (rhs,rhs,n);
00196 fprintf (stdout,"\n norm of residual vector %le",norres);
00197 if (norres>err){
00198
00199
00200 for (j=0;j<ini;j++){
00201
00202
00203 Mp->ssle->solve_system (Gtm,Smat,lhs,rhs,Out);
00204
00205 for (k=0;k<n;k++){
00206 lhs[k]+=p[k];
00207 a[k] = 1.0/dt/dt/beta*(lhs[k]-dd[k]);
00208 v[k] = vv[j] + gamma*dt*a[k];
00209 }
00210
00211
00212 copyv (brhs,rhs,n);
00213
00214
00215
00216 Mmat->gmxv (a,p);
00217 cmulv (-1.0*dt*dt*beta,p,n);
00218 addv (rhs,p,n);
00219
00220
00221
00222 Dmat->gmxv (v,p);
00223 cmulv (-1.0*dt*dt*beta,p,n);
00224 addv (rhs,p,n);
00225
00226
00227 internal_forces (lcid,p);
00228 cmulv (-1.0*dt*dt*beta,p,n);
00229 addv (rhs,p,n);
00230
00231
00232 norres = ss (rhs,rhs,n);
00233 fprintf (stdout,"\n norm of residual vector %le",norres);
00234 if (norres<err){
00235 break;
00236 }
00237 }
00238
00239 }
00240
00241 if (j==ini){
00242
00243
00244
00245 Mp->time = Mp->timecon.actualtime ();
00246
00247 dt = Mp->timecon.actualforwtimeincr ();
00248 }
00249 else{
00250 compute_req_val (lcid);
00251 print_step(lcid, i, Mp->time, rhs);
00252 print_flush();
00253
00254
00255 Mp->time = Mp->timecon.newtime ();
00256 dt = Mp->timecon.actualbacktimeincr ();
00257
00258
00259 i++;
00260 }
00261
00262 }while(Mp->time<end_time);
00263
00264 print_close();
00265 delete [] p;
00266 delete [] dd;
00267 delete [] vv;
00268 }
00269
00270