00001 #include "plssolver.h"
00002 #include "pglobal.h"
00003 #include "seqfile.h"
00004 #include "genfile.h"
00005 #include <string.h>
00006 #include "mpi.h"
00007
00008 void par_solve_nonlinear_statics ()
00009 {
00010 long i,j,k,n,ni,ini,stop,modif,li;
00011 double a0,a1,a2,d,l1,l2,dl,dlmax,dlmin,psi,lambda,blambda,dlambda,ddlambda,norf,norfp,norv,zero,ierr;
00012 double lambdao,ss1,ss2,ss3,ss4,ss5;
00013 double *r,*ra,*ddr,*u,*v,*f,*fa,*fp,*fc,*fi;
00014
00015 FILE *gr;
00016 gr = fopen (Mp->auxfile,"wt");
00017
00018
00019 n = Ndof;
00020
00021 ni = Mp->nial;
00022
00023 ini = Mp->niilal;
00024
00025 zero = Mp->zero;
00026
00027 ierr = Mp->erral;
00028
00029 dl = Mp->dlal;
00030
00031 dlmax = Mp->dlmaxal;
00032
00033 dlmin = Mp->dlminal;
00034
00035 psi = Mp->psial;
00036
00037
00038 r = new double [n];
00039 ddr = new double [n];
00040 u = new double [n];
00041 v = new double [n];
00042 f = new double [n];
00043 fa = new double [n];
00044 fi = new double [n];
00045
00046
00047 Mm->alloc_backup ();
00048
00049
00050 ra = Lsrs->give_lhs (lcid);
00051 fp = Lsrs->give_rhs (lcid*2);
00052 fc = Lsrs->give_rhs (lcid*2+1);
00053
00054
00055 seldofinit ();
00056
00057 if (Mp->hdbackupal==2){
00058 arclopen (li,n,lambda,dl,ra,fp);
00059 }
00060 else{
00061 lambda=0.0;
00062 lambdao=0.0;
00063 li=0;
00064 }
00065
00066
00067
00068
00069
00070 modif=0;
00071
00072
00073
00074
00075 for (i=li;i<ni;i++){
00076
00077
00078 for (j=0;j<n;j++){
00079 r[j]=ra[j];
00080 }
00081
00082 blambda=lambda;
00083
00084 Mm->backup_state_var ();
00085
00086
00087 aux_nonlin_print (gr,ra,lambda);
00088
00089
00090
00091
00092 fprintf (stdout,"\n arc-length: increment %ld lambda %e dl %e",i,lambda,dl);
00093 if (Mespr==1){
00094
00095
00096
00097
00098
00099
00100
00101
00102
00103
00104 if (Mp->adp==1)
00105 Mp->ado->print (Out,i,lambda,0);
00106
00107 }
00108
00109
00110 stiffness_matrix (lcid);
00111
00112
00113 for (j=0;j<n;j++){
00114 f[j]=fp[j];
00115 }
00116
00117
00118 Psol->par_linear_solver (Gt,Smat,v,f,Out,Mespr);
00119
00120
00121
00122
00123
00124
00125
00126
00127
00128
00129
00130
00131
00132
00133
00134 dlambda = dl/sqrt(norv+psi*psi*norfp);
00135
00136
00137
00138
00139 for (j=0;j<n;j++){
00140 ddr[j]=dlambda*v[j];
00141 ra[j]+=ddr[j];
00142 fa[j]=fc[j]+(lambda+dlambda)*fp[j];
00143 }
00144
00145
00146
00147 ddlambda=dlambda;
00148
00149
00150 internal_forces (lcid,fi);
00151
00152 for (k=0;k<n;k++){
00153 f[k]=fa[k]-fi[k];
00154 }
00155
00156 norf=ss(f,f,n);
00157 norf=sqrt(norf);
00158
00159 if (Mespr==1) fprintf (stdout,"\n %e %e norf %e",lambda,dl,norf);
00160
00161 if (norf<ierr){
00162
00163
00164
00165 modif++;
00166
00167 if (modif>1){
00168
00169 dl*=2.0;
00170 if (dl>dlmax) dl=dlmax;
00171 modif=0;
00172
00173 if (Mespr==1){
00174 fprintf (stdout,"\n arc length must be modified (dl increase) dl=%e norf=%e",dl,norf);
00175 fprintf (Out,"\n arc length must be modified (dl increase) dl=%e norf=%e",dl,norf);
00176 }
00177 }
00178
00179 lambda+=dlambda;
00180 continue;
00181 }
00182
00183
00184
00185
00186 stop=0;
00187 for (j=0;j<ini;j++){
00188
00189
00190
00191 Psol->par_linear_solver (Gt,Smat,u,f,Out,Mespr);
00192
00193
00194 for (k=0;k<n;k++){
00195 f[k]=ddr[k];
00196 ddr[k]+=u[k];
00197 }
00198
00199
00200 a2=norv+psi*psi*norfp;
00201 a1=2.0*ss(ddr,v,n)+2.0*ddlambda*psi*psi*norfp;
00202 a0=ss(ddr,ddr,n)+ddlambda*ddlambda*psi*psi*norfp-dl*dl;
00203
00204
00205
00206
00207
00208
00209
00210
00211
00212
00213
00214
00215
00216
00217
00218
00219
00220
00221 if (fabs(a2)<zero){
00222 if (fabs(a1)<zero){
00223 if (fabs(a0)>zero) fprintf (stderr,"\n\n nonsolvable constrained condition in function arclength");
00224 else fprintf (stderr,"\n\n infinite number of solution of constrained condition in function arclength");
00225 }
00226 else{
00227 dlambda=(0.0-a0)/a1;
00228 }
00229 }
00230 else{
00231 d=a1*a1-4.0*a2*a0;
00232
00233 if (d<0.0){
00234 fprintf (stderr,"\n negative discriminant in function arclength");
00235 break;
00236 }
00237
00238
00239 l1=(0.0-a1+sqrt(d))/2.0/a2;
00240 l2=(0.0-a1-sqrt(d))/2.0/a2;
00241
00242
00243
00244 }
00245
00246
00247
00248 ss1=0.0; ss2=0.0;
00249 ss3=0.0; ss4=0.0; ss5=0.0;
00250 for (k=0;k<n;k++){
00251 ss1+=(ddr[k]+l1*v[k])*f[k];
00252 ss2+=(ddr[k]+l2*v[k])*f[k];
00253 ss3+=(ddr[k]+l1*v[k])*(ddr[k]+l1*v[k]);
00254 ss4+=(ddr[k]+l2*v[k])*(ddr[k]+l2*v[k]);
00255 ss5+=f[k]*f[k];
00256 }
00257
00258 if (ss1/ss3/ss5>ss2/ss4/ss5) dlambda=l1;
00259 else dlambda=l2;
00260
00261
00262 for (k=0;k<n;k++){
00263 ddr[k]+=dlambda*v[k];
00264 ra[k]+=u[k]+dlambda*v[k];
00265 fa[k]+=dlambda*fp[k];
00266 }
00267 ddlambda+=dlambda;
00268
00269
00270 internal_forces (lcid,fi);
00271
00272 for (k=0;k<n;k++){
00273 f[k]=fa[k]-fi[k];
00274 }
00275
00276 norf=ss(f,f,n);
00277 norf=sqrt(norf);
00278
00279 if (Mespr==1){
00280 fprintf (stdout,"\n increment %ld inner loop %ld norf=%e",i,j,norf);
00281 fprintf (Out,"\n increment %ld inner loop %ld norf=%e",i,j,norf);
00282 }
00283
00284 if (norf<ierr){
00285 lambda+=ddlambda;
00286 stop=1; break;
00287 }
00288 }
00289
00290 modif=0;
00291
00292 if (stop==0){
00293
00294 dl/=2.0;
00295 if (dl<dlmin){ dl=dlmin; break; }
00296
00297 if (Mespr==1){
00298 fprintf (stdout,"\n arc length must be modified (dl decrease) dl=%e norf=%e",dl,norf);
00299 fprintf (Out,"\n arc length must be modified (dl decrease) dl=%e norf=%e",dl,norf);
00300 }
00301
00302
00303
00304 for (j=0;j<n;j++){
00305 ra[j]=r[j];
00306 }
00307
00308
00309 lambda=blambda;
00310
00311 Mm->restore_state_var ();
00312
00313
00314 }
00315 }
00316
00317
00318 if (Mp->hdbackupal==1){
00319 arclsave (i,n,lambda,dl,ra,fp);
00320 }
00321
00322
00323 delete [] r;
00324 delete [] fi;
00325 delete [] fa;
00326 delete [] f;
00327 delete [] v;
00328 delete [] u;
00329 delete [] ddr;
00330
00331 fclose (gr);
00332
00333 }
00334
00335
00336
00337
00338
00339
00340 void parseldofinit ()
00341 {
00342 long i,j,k,l,ndofn;
00343 double x1,x2,y1,y2,z1,z2,length;
00344
00345 switch (Mp->nlman->displnorm){
00346 case alldofs:{ break; }
00347 case seldofs:
00348 case seldofscoord:{
00349 for (i=0;i<Mp->nlman->nsdofal;i++){
00350 j=Mp->nlman->seldofal[i];
00351 Mp->nlman->seldofal[i]=Mt->give_dof (Mp->nlman->selnodal[i],j)-1;
00352 }
00353 break;
00354 }
00355 case selnodes:{
00356 Mp->nlman->nsdofal=0;
00357 for (i=0;i<Mp->nlman->nsnal;i++){
00358 Mp->nlman->nsdofal+=Mt->give_ndofn (Mp->nlman->selnodal[i]);
00359 }
00360 Mp->nlman->seldofal = new long [Mp->nlman->nsdofal];
00361 k=0;
00362 for (i=0;i<Mp->nlman->nsnal;i++){
00363 ndofn=Mt->give_ndofn (Mp->nlman->selnodal[i]);
00364 for (j=0;j<ndofn;j++){
00365 l=Mt->give_dof (Mp->nlman->selnodal[i],j);
00366 if (l>0){
00367 Mp->nlman->seldofal[k]=l-1;
00368 k++;
00369 }
00370 }
00371 }
00372 Mp->nlman->nsdofal=k;
00373 break;
00374 }
00375 default:{
00376 fprintf (stderr,"\n\n unknown displacement norm is required in function seldofinit (%s, line %d)",__FILE__,__LINE__);
00377 }
00378 }
00379
00380 }