00001 #include "cpnnpsolvert.h"
00002 #include "globalt.h"
00003 #include "globmatt.h"
00004 #include "elemswitcht.h"
00005 #include "transprint.h"
00006 #include <string.h>
00007
00008
00009
00010
00011
00012
00013
00014 void solve_nonstationary_growing_problem_nonlin ()
00015 {
00016 long i,j,k,nt,tncd,tnce,ini;
00017 double zero,dt,prev_timet,end_time,alpha;
00018 double *d,*p,*lhst,*tdlhst,*rhst;
00019 double *fbt,*fit,*lhst_last;
00020
00021 double norf_last,err,norf;
00022
00023
00024 ini = Tp->nii;
00025
00026 err = Tp->err;
00027
00028 nt=Ndoft;
00029
00030
00031 lhst = Lsrst->give_lhs (0);
00032 tdlhst = Lsrst->give_tdlhs (0);
00033
00034 rhst = Lsrst->give_rhs (0);
00035
00036
00037 d = new double [nt];
00038
00039 p = new double [nt];
00040 fbt = new double [nt];
00041 fit = new double [nt];
00042 lhst_last = new double [nt];
00043
00044
00045 nullv (lhst,nt);
00046 nullv (tdlhst,nt);
00047 nullv (d,nt);
00048 nullv (p,nt);
00049 nullv (fbt,nt);
00050 nullv (fit,nt);
00051
00052 alpha=Tp->alpha;
00053 zero=Tp->zero;
00054
00055
00056 Tp->time=Tp->timecont.starttime ();
00057
00058 dt=Tp->timecont.initialtimeincr ();
00059
00060 end_time = Tp->timecont.endtime ();
00061
00062 prev_timet=Tp->time;
00063
00064
00065 Tm->initmaterialmodels();
00066
00067 approximation ();
00068
00069
00070
00071
00072 i=0;
00073 print_initt(-1, "wt");
00074 print_stept(0,i,Tp->time,NULL);
00075 print_flusht();
00076
00077
00078 Tt->lhs_save (lhst,Lsrst->lhsi,tdlhst);
00079
00080 do{
00081
00082 if (Mesprt != 0) fprintf (stdout,"\n\n increment number %ld, time %f, time step = %f\n",i,Tp->time,dt);
00083
00084
00085
00086
00087 tnce = Gtt->search_newelem (Tp->time,prev_timet);
00088
00089
00090 Gtt->update_elem (Tp->time);
00091
00092
00093
00094
00095 tncd = Gtt->search_newdofs (Tp->time,prev_timet);
00096
00097
00098 Gtt->update_dofs (Tp->time);
00099
00100
00101 Ndoft = Gtt->codenum_generation (Outt);
00102 nt=Ndoft;
00103
00104 Tm->updateipval ();
00105 i++;
00106
00107 if (tnce!=0 || tncd!=0){
00108
00109 Tt->lhs_restore (lhst,Lsrst->lhsi,tdlhst);
00110 }
00111
00112
00113 if (Kmat != NULL){
00114 delete Kmat;
00115 Kmat=NULL;
00116 }
00117
00118 if (Cmat != NULL){
00119 delete Cmat;
00120 Cmat=NULL;
00121 }
00122
00123
00124
00125
00126
00127 if (tnce>0){
00128 Tt->initial_nodval ();
00129 }
00130
00131
00132 conductivity_matrix (0);
00133
00134
00135 capacity_matrix (0);
00136
00137 for (j=0;j<nt;j++){
00138 p[j] = lhst[j] + (1.0-alpha)*dt*tdlhst[j];
00139 }
00140
00141
00142 trfel_right_hand_side (0,rhst,nt);
00143
00144
00145 Kmat->gmxv (p,d);
00146
00147
00148
00149 Kmat->scalgm (dt*alpha);
00150 Kmat->addgm (1.0,*Cmat);
00151
00152 for (j=0;j<nt;j++){
00153 rhst[j] = rhst[j] - d[j];
00154 }
00155
00156
00157 Tp->ssle->solve_system (Gtt,Kmat,tdlhst,rhst,Outt);
00158
00159 for (j=0;j<nt;j++){
00160
00161
00162
00163
00164 lhst[j] = p[j] + alpha*dt*tdlhst[j];
00165 }
00166
00167
00168 solution_correction ();
00169
00170 Tt->lhs_save (lhst,Lsrst->lhsi,tdlhst);
00171
00172
00173 approximation ();
00174
00175 norf_last = 1.0e20;
00176
00177 for (j=0;j<ini;j++){
00178
00179
00180 solution_correction ();
00181
00182 Tt->lhs_save (lhst,Lsrst->lhsi,tdlhst);
00183
00184
00185 approximation ();
00186
00187
00188 if (Tp->trsolv == fullnewtont){
00189
00190
00191 capacity_matrix (0);
00192
00193
00194 conductivity_matrix (0);
00195
00196
00197 Kmat->gmxv (p,d);
00198
00199
00200
00201 Kmat->scalgm (dt*alpha);
00202 Kmat->addgm (1.0,*Cmat);
00203 }
00204
00205 if (Tp->trestype==lrhst){
00206
00207 trfel_right_hand_side (0,rhst,nt);
00208
00209 Kmat->gmxv (tdlhst,fit);
00210
00211 for (k=0;k<nt;k++){
00212 fbt[k] = rhst[k] - d[k] - fit[k];
00213 }
00214 }
00215
00216 if (Tp->trestype==fluxest){
00217
00218 internal_fluxes (fit,nt);
00219
00220 for (k=0;k<nt;k++){
00221 fbt[k]=fit[k];
00222 }
00223 }
00224
00225 norf = ss (fbt,fbt,nt);
00226
00227 if (Mesprt==1) fprintf (stdout,"\n inner iteration %ld error %e",j,norf);
00228
00229 if (norf<err) break;
00230
00231 Tp->ssle->solve_system (Gtt,Kmat,d,fbt,Outt);
00232
00233 for (k=0;k<nt;k++){
00234 tdlhst[k]+=d[k];
00235 lhst[k]+=alpha*dt*d[k];
00236 }
00237
00238
00239 solution_correction ();
00240
00241 Tt->lhs_save (lhst,Lsrst->lhsi,tdlhst);
00242
00243
00244 approximation ();
00245
00246
00247
00248 if (Tp->convergcontrolt==yes){
00249 if (norf > norf_last){
00250 for (k=0;k<nt;k++){
00251 lhst[k] = lhst_last[k];
00252 }
00253
00254
00255 solution_correction ();
00256
00257 Tt->lhs_save (lhst,Lsrst->lhsi,tdlhst);
00258
00259
00260 approximation ();
00261
00262 if (Mesprt==1) fprintf (stdout,"\n\n convergence control: inner iteration skiped %ld error %e\n\n",j,norf);
00263
00264 break;
00265 }
00266 for (k=0;k<nt;k++){
00267 lhst_last[k] = lhst[k];
00268 }
00269
00270 norf_last = norf;
00271 }
00272 }
00273
00274
00275 prev_timet=Tp->time;
00276 Tp->time = Tp->timecont.newtime ();
00277 dt = Tp->timecont.actualbacktimeincr ();
00278
00279 print_stept(0,i,Tp->time,NULL);
00280 print_flusht();
00281
00282 }while(Tp->time<=end_time);
00283
00284 print_closet ();
00285
00286 delete [] p;
00287 delete [] d;
00288
00289 delete [] fit;
00290 delete [] fbt;
00291 delete [] lhst_last;
00292 }