00001 #include "dnpsolvert.h"
00002 #include "globalt.h"
00003 #include "globmatt.h"
00004 #include "transprint.h"
00005 #include "saddpoint.h"
00006 #include <string.h>
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016 void solve_discont_nonstationary_problem ()
00017 {
00018 lin_nonstat_dform ();
00019
00020
00021
00022 }
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032 void lin_nonstat_dform ()
00033 {
00034 long i,j,n,nbdof,nsad,lcid;
00035 double dt,end_time,alpha,zero,totsalt;
00036 double *f,*d,*p,*lhs,*tdlhs,*rhs;
00037
00038
00039 lcid=0;
00040
00041
00042
00043 n=Ndoft;
00044
00045 nbdof=Gtt->nbdof;
00046
00047
00048 nsad=Gtt->nsad;
00049
00050
00051 zero=Tp->zero;
00052
00053
00054
00055 lhs = Lsrst->give_lhs (0);
00056 nullv (lhs,n);
00057
00058 tdlhs = Lsrst->give_tdlhs (0);
00059 nullv (tdlhs,n);
00060
00061 rhs = Lsrst->give_rhs (0);
00062 nullv (rhs,n);
00063
00064
00065 f = new double [n];
00066 nullv (f,n);
00067
00068 d = new double [n];
00069 nullv (d,n);
00070
00071 p = new double [n];
00072 nullv (p,n);
00073
00074
00075 alpha=Tp->alpha;
00076
00077
00078 Tp->time=Tp->timecont.starttime ();
00079
00080 end_time = Tp->timecont.endtime ();
00081
00082 dt = Tp->timecont.initialtimeincr ();
00083
00084 i=0;
00085
00086
00087 approximation ();
00088
00089
00090 if (Tp->nvs==1 && Tp->pnvs==1)
00091 actual_previous_nodval ();
00092
00093
00094
00095 Tm->initmaterialmodels();
00096
00097 print_initt(-1, "wt");
00098 print_stept(0,i,Tp->time,NULL);
00099
00100
00101
00102
00103 do{
00104
00105
00106 Tp->time=Tp->timecont.newtime (dt);
00107
00108 if (Mesprt != 0) fprintf (stdout,"\n\n iteration number %ld, time %f, time step = %f\n",i,Tp->time,dt);
00109
00110 Tm->updateipval ();
00111
00112 i++;
00113
00114
00115 capacity_matrix (0);
00116
00117
00118 conductivity_matrix (0);
00119
00120
00121
00122 for (j=0;j<n;j++){
00123 d[j] = lhs[j]+dt*(1.0-alpha)*tdlhs[j];
00124 }
00125
00126
00127 Cmat->gmxv (d,p);
00128
00129
00130
00131 Kmat->scalgm (dt*alpha);
00132 Kmat->addgm (1.0,*Cmat);
00133
00134
00135
00136 trfel_right_hand_side (lcid,f,n);
00137
00138
00139 for (j=0;j<n;j++){
00140 rhs[j] = f[j]*alpha*dt + p[j];
00141 }
00142
00143 Tt->compute_jumps (rhs);
00144
00145
00146 Gtt->mult_localization (Kmat);
00147
00148
00149 Kmat->diag_check (zero,rhs);
00150
00151
00152
00153
00154 Tp->ssle->solve_system (Gtt,Kmat,lhs,rhs,Outt);
00155
00156
00157
00158
00159 for (j=0;j<n;j++){
00160 tdlhs[j]=(lhs[j]-d[j])/dt/alpha;
00161 }
00162
00163
00164 solution_correction ();
00165
00166 approximation ();
00167
00168 totsalt = total_integral_ip (5);
00169 fprintf (stdout,"\n total salt content %lf",totsalt);
00170 totsalt = total_integral_ip (6);
00171 fprintf (stdout,"\n integral of cf %lf",totsalt);
00172 totsalt = total_integral_ip (7);
00173 fprintf (stdout,"\n integral of cc %lf",totsalt);
00174 totsalt = total_integral_ip (8);
00175 fprintf (stdout,"\n integral of cb %lf",totsalt);
00176
00177 if (Tp->nvs==1 && Tp->pnvs==1)
00178 actual_previous_nodval ();
00179
00180
00181 print_stept(0,i,Tp->time,NULL);
00182 print_flusht();
00183
00184 }while(Tp->time<end_time);
00185
00186 delete [] p;
00187 delete [] d;
00188 delete [] f;
00189
00190 print_closet();
00191 }
00192
00193
00194
00195
00196
00197
00198
00199
00200
00201
00202
00203
00204
00205
00206
00207
00208
00209
00210
00211
00212
00213
00214
00215
00216 void lin_nonstat_dform_resistance ()
00217 {
00218 long i,j,n,ii,ipp,nbdof,nsad,lcid;
00219 double dt,end_time,alpha,zero;
00220 double *f,*d,*p,*lhs,*tdlhs,*rhs;
00221
00222
00223 lcid=0;
00224
00225
00226
00227 n=Ndoft;
00228
00229 nbdof=Gtt->nbdof;
00230
00231
00232 nsad=Gtt->nsad;
00233
00234
00235 zero=Tp->zero;
00236
00237
00238
00239 lhs = Lsrst->give_lhs (0);
00240 nullv (lhs,n);
00241
00242 tdlhs = Lsrst->give_tdlhs (0);
00243 nullv (tdlhs,n);
00244
00245 rhs = Lsrst->give_rhs (0);
00246 nullv (rhs,n);
00247
00248
00249 f = new double [n];
00250 nullv (f,n);
00251
00252 d = new double [n];
00253 nullv (d,n);
00254
00255 p = new double [n];
00256 nullv (p,n);
00257
00258
00259 alpha=Tp->alpha;
00260
00261
00262 Tp->time=Tp->timecont.starttime ();
00263
00264 end_time = Tp->timecont.endtime ();
00265
00266 dt = Tp->timecont.initialtimeincr ();
00267
00268 i=0;
00269
00270
00271 approximation ();
00272
00273
00274 if (Tp->nvs==1 && Tp->pnvs==1)
00275 actual_previous_nodval ();
00276
00277
00278
00279 Tm->initmaterialmodels();
00280
00281 print_initt(-1, "wt");
00282 print_stept(0,i,Tp->time,NULL);
00283
00284
00285
00286
00287 do{
00288
00289
00290 Tp->time=Tp->timecont.newtime (dt);
00291
00292 if (Mesprt != 0) fprintf (stdout,"\n\n iteration number %ld, time %f, time step = %f\n",i,Tp->time,dt);
00293
00294 Tm->updateipval ();
00295
00296 i++;
00297
00298
00299 capacity_matrix (0);
00300
00301
00302 conductivity_matrix (0);
00303
00304
00305
00306 for (j=0;j<n;j++){
00307 d[j] = lhs[j]+dt*(1.0-alpha)*tdlhs[j];
00308 }
00309
00310
00311 Cmat->gmxv (d,p);
00312
00313
00314
00315 Kmat->scalgm (dt*alpha);
00316 Kmat->addgm (1.0,*Cmat);
00317
00318
00319
00320 trfel_right_hand_side (lcid,f,n);
00321
00322
00323 for (j=0;j<n;j++){
00324 rhs[j] = f[j]*alpha*dt + p[j];
00325 }
00326
00327 Tt->compute_jumps (rhs);
00328
00329
00330 Gtt->mult_localization (Kmat);
00331
00332
00333
00334
00335
00336 Tp->ssle->solve_system (Gtt,Kmat,lhs,rhs,Outt);
00337
00338
00339
00340
00341 for (j=0;j<n;j++){
00342 tdlhs[j]=(lhs[j]-d[j])/dt/alpha;
00343 }
00344
00345
00346 solution_correction ();
00347
00348 approximation ();
00349
00350
00351 for (ii=0;i<Tt->ne;ii++){
00352
00353 ipp=Tt->elements[ii].ipp[0][0];
00354
00355 Tm->computenlfluxes (0,ipp);
00356 Tm->computenlfluxes (1,ipp);
00357 }
00358 Tt->compute_resistance_factor (rhs);
00359
00360
00361
00362
00363
00364 if (Tp->nvs==1 && Tp->pnvs==1)
00365 actual_previous_nodval ();
00366
00367
00368 print_stept(0,i,Tp->time,NULL);
00369 print_flusht();
00370
00371 }while(Tp->time<end_time);
00372
00373 delete [] p;
00374 delete [] d;
00375 delete [] f;
00376
00377 print_closet();
00378 }
00379