00001 #include "arclength.h" 00002 #include "nssolver.h" 00003 #include "mechprint.h" 00004 #include "mathem.h" 00005 #include "global.h" 00006 #include "globmat.h" 00007 #include "loadcase.h" 00008 00009 #include "mtsolver.h" 00010 00011 #include <stdlib.h> 00012 00013 00014 00015 /** 00016 The function solves system of nonlinear algebraic equations 00017 by the arc-length method. This method was modified to reach the given value lambdar 00018 of the lambda parameter. lambdar value is reached with required error errl. 00019 Only one right hand side %vector is supported with respect to nonlinearity and absence of superposition. 00020 Proportional %vector grows until the nial number of steps is reached. 00021 00022 @param lcid - load case id 00023 @param nlman - pointer to structure conatining setup of the solver 00024 @param ra - %vector of attained displacements 00025 @param fa - attained load %vector 00026 @param fc - constant part of load %vector (it is not influenced by lambda) 00027 @param fp - proportional part of load %vector (it is increased by load coefficient lambda up to the nlman->lambdar) 00028 @param flc - constant component of load %vector due to forces 00029 @param flp - proportional component of load %vector due to forces 00030 @param li - initial value of step id (default is 0) 00031 @param ilambda - initial value of load coefficient (default is 0) 00032 @param outres - flag for performing of output of results (if yes -> print_step procedure is called for the each step) 00033 00034 @return The function returns reached lambda parameter. 00035 00036 Created by JK, 16.8.2001 00037 Rewritten by TKo, JK 08.2011 00038 */ 00039 double garclength(long lcid, nonlinman *nlman, double *ra, double *fa, double *fc, double *fp, double *flc, double *flp, 00040 long li, double ilambda, answertype outres) 00041 { 00042 long i,j,n,ni,ini,stop,stopj,modif; 00043 double dl; // actual value of arc length increment 00044 double dlmin; // minimum value of arc length increment 00045 double dlmax; // maximum value of arc length increment 00046 double psi; // proportional coefficient between load vector and displacement vector 00047 double psi0; // initial value of proportional coefficient psi 00048 double lambda; // actual value of load coefficient 00049 double blambda; // backupped value of lambda 00050 double dlambda; // load coefficient increment from the inner loop 00051 double ddlambda; // cummulative value of load coefficent increments ddlambda = \sum(dlambda) 00052 double ddlambda0;// initial/reference value of ddlambda 00053 double norf; // norm of residual vector 00054 double norfp; // norm of proportional load vector 00055 double norv; // norm of displacement increments 00056 double norfa; // norm of attained load vector 00057 double zero; // computer zero 00058 double ierr; // relative error of residual vector 00059 double dtr; // time step size reduction coeffcient required by material models 00060 double tmp; // auxiliary variable 00061 double *r; // displacements from previous step (backup copy) 00062 double *ddr; // vector of displacement increments 00063 double *u; // vector of corrections for computed displacements v 00064 double *v; // vector of computed displacements for actual arclength step 00065 double *f; // auxiliary vector of actual load (call of solve_system rewrite it during solution) 00066 double *fi; // vector of internal forces 00067 double totl = 0.0; // the total(cumulative) length of arc 00068 long sm_change = 0; // flag for changed content of the stiffness matrix 00069 long modif_psi; // flag for automatic modification psi in dependence on the ddlambda 00070 flagsw check_rv_fp = nlman->check_lambdar; // flag for checking of attainment of required value for load coefficient lambdar 00071 flagsw check_max_tot_al = nlman->check_tot_al; // flag fro checking maximu value of the length of arc 00072 matrix lsm_a(3,3); // least square matrix for the divergence detection 00073 vector lsm_r(3), lsm_l(3); // right hand and left hand side for the divergence detection 00074 00075 00076 00077 // number of rows of the matrix 00078 n = Ndofm; 00079 // computer zero 00080 zero = Mp->zero; 00081 // initial value of load coefficient 00082 lambda = ilambda; 00083 // maximum number of increments 00084 ni = nlman->nial; 00085 // maximum number of iterations in one increment 00086 ini = nlman->niilal; 00087 // required error in inner loop 00088 ierr = nlman->erral; 00089 // length of arc 00090 dl = nlman->dlal; 00091 // maximum length of arc 00092 dlmax = nlman->dlmaxal; 00093 // minimum length of arc 00094 dlmin = nlman->dlminal; 00095 // displacement-loading driving switch 00096 psi = nlman->psial; 00097 00098 // allocation of auxiliary arrays 00099 r = new double[n]; 00100 ddr = new double[n]; 00101 u = new double[n]; 00102 v = new double[n]; 00103 f = new double[n]; 00104 fi = new double[n]; 00105 00106 // Delta lambda - cumulative increment in one load increment 00107 // at the beginning of computation is always zero 00108 dlambda=0.0; 00109 00110 /* 00111 // assemble initial load vector, fa[j]=fc[j]+lambda*fp[j] 00112 // in case that no restorage from the backup was performed (lambda == 0.0) 00113 if (lambda == 0.0) 00114 assemble_attained_load_vector(fa, fc, fp, n, lambda, nlman); 00115 */ 00116 00117 // in case that no restorage from the backup was performed (lambda == 0.0) 00118 if (lambda == 0.0) 00119 { 00120 sm_change = assemble_stiffness_matrix (lcid,li,-1,li); 00121 // backup of the fc because the right hand side will be destroyed during the solution of the system of equations 00122 copyv(fc, f, n); 00123 // solution of K(r).v=F 00124 Mp->ssle->solve_system(Gtm,Smat,ra,f,Out); 00125 } 00126 00127 // norm of proportionality vector due to forces 00128 // norfp = normv(fp, n); 00129 norfp = normv(flp, n); 00130 // norfp = loadincr(fp,n); 00131 modif=0; 00132 psi0 = psi; 00133 if (psi<0.0){ 00134 // parameter psi will be modified during the computation 00135 modif_psi=1; 00136 }else{ 00137 // parameter psi will be constant during the computation 00138 modif_psi=0; 00139 } 00140 00141 00142 // *************************** 00143 // main iteration loop **** 00144 // *************************** 00145 for (i=li;i<ni;i++) 00146 { // loop i 00147 00148 // indicator of strain computation 00149 // it means, strains at integration points have not been computed 00150 Mp->strainstate=0; 00151 // indicator of stress computation 00152 // it means, stresses at integration points have not been computed 00153 Mp->stressstate=0; 00154 // indicator of computation of other array 00155 // it means, stresses at integration points have not been computed 00156 Mp->otherstate=0; 00157 00158 Mp->istep = i; 00159 Mp->jstep = -1; 00160 Mp->lambda = lambda; 00161 Mp->dlambda = 0.0; 00162 00163 // backup of left hand side vector rb[j]=ra[j]; 00164 copyv(ra, r, n); 00165 // backup of reached lambda parameter 00166 blambda=lambda; 00167 00168 if (Mespr==1) 00169 { 00170 fprintf (stdout,"\n increment %ld, lambda=%e dl=%e psi=%e",i,lambda,dl,psi); 00171 // fprintf (Out,"\n increment %ld, lambda=%e dl=%e",i,lambda,dl); 00172 } 00173 00174 // assembling of the stiffness matrix 00175 sm_change = assemble_stiffness_matrix (lcid,i,-1,li); 00176 00177 // backup of the fp because the right hand side will be destroyed during the solution of the system of equations 00178 copyv(fp, f, n); 00179 00180 // solution of K(r).v=F 00181 Mp->ssle->solve_system(Gtm,Smat,v,f,Out); 00182 00183 // square of generalized norm of displacement increments 00184 norv = displincr(v,v,n,lcid); 00185 00186 if ((i == 0) && (psi < 0.0) && (norfp > Mp->zero)) 00187 psi0 = psi = sqrt(norv/norfp/norfp); 00188 00189 // compute new dlambda increment 00190 tmp = norv+psi*psi*norfp*norfp; 00191 // compute new dlambda increment 00192 dlambda = dl/sqrt(tmp); 00193 00194 // check required value of load coefficient lambdar 00195 if ((check_rv_fp == on) && (dlambda+lambda > nlman->lambdar)) 00196 { 00197 check_rv_fp = off; 00198 dlambda = nlman->lambdar - lambda; 00199 dl = dlambda*sqrt(tmp); 00200 } 00201 // checking of maximum value of the total (cumulative) arc length 00202 if ((check_max_tot_al == on) && (dl+totl > nlman->max_tot_al)) 00203 { 00204 check_max_tot_al = off; 00205 dl = nlman->max_tot_al - totl; 00206 } 00207 00208 // ddr[j]=dlambda*v[j]; 00209 copymultv(v, ddr, dlambda, n); 00210 00211 // ra[j]+=ddr[j]; 00212 addv(ra, ddr, n); 00213 00214 /* 00215 // update attained load vector 00216 // fa[j] = fc + (lambda+dlambda)*fp[j] is calculated 00217 addmultv(fc, fp, lambda+dlambda, fa, n); 00218 */ 00219 // update attained load vector 00220 // fa[j] = flc + (lambda+dlambda)*flp[j] is calculated 00221 addmultv(flc, flp, lambda+dlambda, fa, n); 00222 00223 // initialize value of ddlambda before inner iteration loop 00224 ddlambda=dlambda; 00225 00226 if (Mespr == 1) 00227 fprintf(Out, "Step %ld: ddlambda = %e, norv = %e, norfp = %e\n", i, ddlambda, norv, norfp); 00228 00229 Mp->dlambda = ddlambda; 00230 Mp->lambda += ddlambda; // due to prescribed displacements 00231 // computation of internal forces 00232 internal_forces (lcid,fi); 00233 Mp->lambda -= ddlambda; // restore state 00234 00235 // f[j] = fa[j]-fi[j] 00236 subv(fa, fi, f, n); 00237 00238 // computation of norms 00239 norf=normv(f,n); 00240 norfa = normv(fa, n); 00241 00242 // proportional norm is used in case of nonzero load vector 00243 // otherwise absolute norm is used for vector of unbalanced forces 00244 if (fabs(norfa) > Mp->zero) 00245 norf /= norfa; 00246 00247 // check time step size requirements which originate in material models, 00248 // after check, the dtr will contain the minimum ratio of required time step size to the actual one 00249 // ( dtr \in (0.0; 1.0) ) collected from all integration points or dtr = 1.0 for no change 00250 // in time step size 00251 dtr = Mm->dstep_red_mat(); 00252 if (dtr < 1.0) // reduction of time step size was required by some material model 00253 { 00254 norf = 2.0*ierr; 00255 if (dl == dlmin) // dl cannot be decreased 00256 break; 00257 // modification of the arc length 00258 dl*=dtr; 00259 if (dl<dlmin) 00260 dl=dlmin; 00261 00262 if (Mespr==1) 00263 { 00264 if (dl == dlmin) 00265 { 00266 fprintf (stdout,"\n arc length was decreased due to mat. models dl=dlmin=%e",dl); 00267 // fprintf (Out,"\n arc length was decreased dl=dlmin=%e",dl); 00268 } 00269 else 00270 { 00271 fprintf (stdout,"\n arc length was decreased due to mat. models dl=%e",dl); 00272 // fprintf (Out,"\n arc length was decreased dl=%e",dl); 00273 } 00274 } 00275 // restoring of left hand side vector 00276 copyv (r, ra, n); 00277 // restoring of lambda parameter 00278 lambda=blambda; 00279 continue; 00280 } 00281 00282 if (norf<ierr) 00283 { 00284 // ****************************************** 00285 // no inner iteration loop is required *** 00286 // ****************************************** 00287 if (Mespr==1) 00288 fprintf(stdout,"\n increment %ld norf=%e ierr=%e ddlambda=%e dlambda=%e",i,norf,ierr, ddlambda, dlambda); 00289 00290 totl += dl; 00291 00292 modif++; 00293 00294 if (modif>1) 00295 { 00296 // arc length modification 00297 dl*=2.0; 00298 if (dl>dlmax) 00299 dl=dlmax; 00300 modif=0; 00301 if (Mespr==1) 00302 fprintf(stdout,"\n arc-length is modified (dl increase) dl=%e",dl); 00303 } 00304 lambda+=dlambda; 00305 00306 // this part modifies parameter psi 00307 if (modif_psi) 00308 adapt_psi_coeff(i, ddlambda, ddlambda0, psi0, psi); 00309 00310 Mm->updateipval(); 00311 if (outres == yes) 00312 { 00313 compute_req_val(lcid); 00314 print_step(lcid, i+1, lambda, fa); 00315 print_flush(); 00316 } 00317 00318 // required value of lambda was attained for the load vector -> end of arclength 00319 if ((nlman->check_lambdar == on) && (check_rv_fp == off)) 00320 break; 00321 00322 // required total(cumulative) arc length has been attained => end of arclength 00323 if ((nlman->check_tot_al == on) && (check_max_tot_al==off)) 00324 break; 00325 00326 continue; 00327 } 00328 00329 // *************************************** 00330 // * inner iteration loop * 00331 // *************************************** 00332 if (Mespr==1) 00333 fprintf (stdout,"\n inner loop: 1 norf=%e ierr=%e ddlambda=%e dlambda=%e",norf,ierr,ddlambda, dlambda); 00334 stop=0; 00335 stopj=0; 00336 fillm(0.0, lsm_a); 00337 fillv(0.0, lsm_r); 00338 for (j=0;j<ini;j++) 00339 { //loop j 00340 // copy of actual step to problem description 00341 // it is used in connection with printing of output values 00342 Mp->jstep = j; 00343 00344 // eventual update of stiffness matrix 00345 Mespr = 0; // switch off the message printing for assembling of stiffness %matrix 00346 sm_change = assemble_stiffness_matrix (lcid,i,j,li); 00347 Mespr = 1; // switch on the message printing for assembling of stiffness %matrix 00348 00349 // solve displacement increments due to vector of unbalanced forces 00350 Mp->ssle->solve_system(Gtm,Smat,u,f,Out); 00351 if (sm_change) 00352 { 00353 // backup of the fp because the right hand side will be destroyed during the solution of the system of equations 00354 copyv(fp, f, n); 00355 Mp->ssle->solve_system(Gtm,Smat,v,f,Out); 00356 } 00357 00358 // backup vector of displacement increments from the previous step 00359 // f[j] = ddr[j] 00360 copyv(ddr, f, n); 00361 00362 // add computed correction of displacement increments due to unbalanced forces 00363 // ddr[j] = ddr[j]+u[j] 00364 addv(ddr, u, n); 00365 00366 // determination of increment of load parameter 00367 determine_dlambda (ddr, u, v, f, n, ddlambda,psi,norfp,dl,dlambda,stopj, nlman, lcid); 00368 if (stopj) 00369 break; 00370 00371 // check required value of load coefficient lambdar 00372 if ((nlman->check_lambdar == on) && (check_rv_fp == off) && 00373 (dlambda+ddlambda+lambda > nlman->lambdar)) 00374 dlambda = nlman->lambdar - ddlambda - lambda; 00375 00376 // actualize vector of displacements increments 00377 // ddr[j] += dlambda*v[j] 00378 addmultv(ddr, v, dlambda, n); 00379 00380 // actualize vector of attained displacements ra 00381 // ra[j] += u[j] + dlambda*v[j] 00382 addv(ra, u, n); 00383 addmultv(ra, v, dlambda, n); 00384 /* 00385 // update attained load vector 00386 // fa[k]+=dlambda*fp[k] is calculated 00387 addmultv(fa, fp, dlambda, n);*/ 00388 00389 // update attained load vector 00390 // fa[k]+=dlambda*flp[k] is calculated 00391 addmultv(fa, flp, dlambda, n); 00392 //update_attained_load_vector(fa, fp, n, lambda+ddlambda, dlambda, nlman); 00393 00394 ddlambda+=dlambda; 00395 00396 Mp->dlambda = ddlambda; 00397 Mp->lambda += ddlambda; 00398 // computation of internal forces 00399 internal_forces (lcid,fi); 00400 Mp->lambda -= ddlambda; 00401 00402 // check time step size requirements which originate in material models, 00403 // after check, the dtr will contain the minimum ratio of required time step size to the actual one 00404 // ( dtr \in (0.0; 1.0) ) collected from all integration points or dtr = 1.0 for no change 00405 // in time step size 00406 dtr = Mm->dstep_red_mat(); 00407 if (dtr < 1.0) // reduction of time step size was required by some material model 00408 { 00409 stopj = 1; 00410 norf = 2.0*ierr; 00411 break; 00412 } 00413 00414 // computation of residual forces f[k] = fa[k] - fi[k] 00415 subv(fa, fi, f, n); 00416 // norm of the vector of residuals 00417 norf=normv(f,n); 00418 // norm of the vector of attained forces 00419 norfa = normv(fa, n); 00420 if (norfa > zero) 00421 { 00422 // if the attained forces are nonzero, the norm of residuals is scaled by the norm of forces 00423 // this scaling is reasonable because nondimensional quantity is obtained 00424 norf /= norfa; 00425 } 00426 00427 if (Mespr==1) 00428 { 00429 fprintf (stdout,"\n increment %ld inner loop %ld norf=%e",i,j,norf); 00430 // fprintf (Out,"\n increment %ld inner loop %ld norf=%e",i,j,norf); 00431 } 00432 00433 // equilibrium was attained 00434 if (norf<ierr) 00435 { 00436 lambda+=ddlambda; 00437 Mm->updateipval(); 00438 compute_req_val (lcid); 00439 print_step(lcid, i+1, lambda, fa); 00440 print_flush(); 00441 totl += dl; 00442 stop=1; 00443 break; 00444 } 00445 00446 // divergence detection with help of least square method 00447 stopj = check_divergency(nlman, lsm_a, lsm_r, lsm_l, j, norf); 00448 if (stopj) 00449 break; 00450 00451 } // end of j 00452 00453 // there is no modification of the length of arc because at least one iteration was performed 00454 modif=0; 00455 00456 // check of the total(cumulative) attained length of arc 00457 if ((nlman->check_tot_al == on) && (check_max_tot_al==off)) 00458 break; 00459 00460 if (stop==1) // the equilibrium was attained 00461 { 00462 // this part modifies parameter psi 00463 if (modif_psi) 00464 adapt_psi_coeff(i, ddlambda, ddlambda0, psi0, psi); 00465 00466 if ((nlman->check_lambdar == on) && (check_rv_fp == off) && (fabs(lambda-nlman->lambdar) > nlman->errl)) 00467 check_rv_fp = on; 00468 00469 } 00470 00471 if ((stop==0) || (stopj)) // stop == 0 the equilibrium was not attained in the prescribed number of iterations 00472 // stopj != 0 the dlambda could not be determined or divergence detected or 00473 // material models required time step reduction 00474 { 00475 if (dl == dlmin) // dl cannot be decreased 00476 break; 00477 // modification of the arc length 00478 if (dtr < 1.0) // reduction of time step size was required by some material model 00479 dl *= dtr; 00480 else 00481 dl/=2.0; 00482 00483 if (dl<dlmin) 00484 dl=dlmin; 00485 00486 if (Mespr==1) 00487 { 00488 if (dl == dlmin) 00489 { 00490 if (dtr < 1.0) 00491 fprintf (stdout,"\n arc length was decreased due to mat. models dl=dlmin=%e",dl); 00492 else 00493 fprintf (stdout,"\n arc length was decreased dl=dlmin=%e",dl); 00494 // fprintf (Out,"\n arc length was decreased dl=dlmin=%e",dl); 00495 } 00496 else 00497 { 00498 if (dtr < 1.0) 00499 fprintf (stdout,"\n arc length was decreased due to mat. models dl=%e",dl); 00500 else 00501 fprintf (stdout,"\n arc length was decreased dl=%e",dl); 00502 // fprintf (Out,"\n arc length was decreased dl=%e",dl); 00503 } 00504 } 00505 // restoring of left hand side vector 00506 copyv (r, ra, n); 00507 // restoring of lambda parameter 00508 lambda=blambda; 00509 } 00510 00511 // ------------------------------------ 00512 // finish of main iteration loop ---- 00513 // ------------------------------------ 00514 } 00515 00516 fprintf(stdout, "\n The total attained lambda =% le\n", lambda); 00517 00518 if ((i < ni-1) && (norf >= ierr)) 00519 { 00520 if (Mespr==1) 00521 { 00522 print_err("the step length was decreased to the minimum required value\n" 00523 " but the equilibrium could not be attained.\n" 00524 " FORCED output of the attained results was performed in this step.", 00525 __FILE__, __LINE__, __func__); 00526 } 00527 compute_req_val (lcid); 00528 print_step_forced(lcid, i+1, lambda, fa); 00529 print_flush(); 00530 } 00531 00532 if (Mp->hdbcont.save_stat()) 00533 { 00534 // backup of the attained state is required 00535 if (Mespr==1) 00536 fprintf (stdout,"\n Creating backup file\n"); 00537 solver_save (ra, fa, i, lambda, dl, NULL, n); 00538 } 00539 00540 00541 delete [] r; 00542 delete [] fi; 00543 delete [] f; 00544 delete [] v; 00545 delete [] u; 00546 delete [] ddr; 00547 00548 return lambda; 00549 } 00550 00551 00552 00553 /** 00554 The function solves system of nonlinear algebraic equations 00555 by the arc-length method. This method was modified to reach the given value lambdar 00556 of the lambda parameter. lambdar value is reached with required error errl. 00557 Only one right hand side %vector is supported with respect to nonlinearity and absence of superposition. 00558 Proportional %vector grows until the nial number of steps is reached. 00559 00560 @param lcid - load case id 00561 @param nlman - pointer to structure conatining setup of the solver 00562 @param ra - %vector of attained displacements 00563 @param fa - attained load %vector 00564 @param fc - constant part of load %vector (it is not influenced by lambda) 00565 @param fp - proportional part of load %vector (it is increased by load coefficient lambda up to the nlman->lambdar) 00566 @param li - initial value of step id (default is 0) 00567 @param ilambda - initial value of load coefficient (default is 0) 00568 @param outres - flag for performing of output of results (if yes -> print_step procedure is called for the each step) 00569 00570 @return The function returns reached lambda parameter. 00571 00572 Created by JK, 16.8.2001 00573 Rewritten by TKo, JK 08.2011 00574 */ 00575 /* 00576 double garclength2(long lcid, nonlinman *nlman, double *ra, double *fa, double *fc, double *fp, 00577 long li, double ilambda, answertype outres) 00578 { 00579 long i,j,n,ni,ini,stop,stopj,modif; 00580 double dl; // actual value of arc length increment 00581 double dlmin; // minimum value of arc length increment 00582 double dlmax; // maximum value of arc length increment 00583 double psi0; // initial value of proportional coefficient psi 00584 double blambda; // backupped value of lambda 00585 double zero; // computer zero 00586 double ierr; // relative error of residual vector 00587 double *r; // displacements from previous step (backup copy) 00588 double *ddr; // vector of displacement increments 00589 double *u; // vector of corrections for computed displacements v 00590 double *v; // vector of computed displacements for actual arclength step 00591 double *f; // auxiliary vector of actual load (call of solve_system rewrite it during solution) 00592 double *fi; // vector of internal forces 00593 double totl = 0.0; // the total(cumulative) length of arc 00594 long sm_change = 0; // flag for changed content of the stiffness matrix 00595 long modif_psi; // flag for automatic modification psi in dependence on the ddlambda 00596 matrix lsm_a(3,3); // least square matrix for the divergence detection 00597 vector lsm_r(3), lsm_l(3); // right hand and left hand side for the divergence detection 00598 flagsw check_rv_fp = nlman->check_lambdar; // flag for checking of attainment of required value for load coefficient lambdar 00599 flagsw check_max_tot_al = nlman->check_tot_al; // flag fro checking maximu value of the length of arc 00600 00601 00602 00603 // number of rows of the matrix 00604 n = Ndofm; 00605 // computer zero 00606 zero = Mp->zero; 00607 // initial value of load coefficient 00608 lambda = ilambda; 00609 // maximum number of increments 00610 ni = nlman->nial; 00611 // displacement-loading driving switch 00612 psi = nlman->psial; 00613 00614 // allocation of auxiliary arrays 00615 r = new double[n]; 00616 ddr = new double[n]; 00617 u = new double[n]; 00618 v = new double[n]; 00619 f = new double[n]; 00620 fi = new double[n]; 00621 00622 00623 00624 // assemble initial load vector, fa[j]=fc[j]+lambda*fp[j] 00625 // in case that no restorage from the backup was performed (lambda == 0.0) 00626 if (lambda == 0.0) 00627 assemble_attained_load_vector(fa, fc, fp, n, lambda, nlman); 00628 // norm of proportionality vector 00629 norfp = normv(fp, n); 00630 // norfp = loadincr(fp,n); 00631 modif=0; 00632 psi0 = psi; 00633 if (psi<0.0){ 00634 // parameter psi will be modified during the computation 00635 modif_psi=1; 00636 }else{ 00637 // parameter psi will be constant during the computation 00638 modif_psi=0; 00639 } 00640 00641 00642 // *************************** 00643 // main iteration loop **** 00644 // *************************** 00645 for (i=li;i<ni;i++) 00646 { // loop i 00647 00648 // indicator of strain computation 00649 // it means, strains at integration points have not been computed 00650 Mp->strainstate=0; 00651 // indicator of stress computation 00652 // it means, stresses at integration points have not been computed 00653 Mp->stressstate=0; 00654 // indicator of computation of other array 00655 // it means, stresses at integration points have not been computed 00656 Mp->otherstate=0; 00657 00658 Mp->istep = i; 00659 Mp->jstep = -1; 00660 00661 // backup of left hand side vector rb[j]=ra[j]; 00662 copyv(ra, r, n); 00663 // backup of reached lambda parameter 00664 blambda=lambda; 00665 00666 if (Mespr==1) 00667 { 00668 fprintf (stdout,"\n increment %ld, lambda=%e dl=%e psi=%e",i,lambda,dl,psi); 00669 // fprintf (Out,"\n increment %ld, lambda=%e dl=%e",i,lambda,dl); 00670 } 00671 00672 // backup of the fp because the right hand side will be destroyed during the solution of the system of equations 00673 copyv(fp, f, n); 00674 00675 00676 // there is no modification of the length of arc because at least one iteration was performed 00677 modif=0; 00678 00679 norf = graclength_one_step() 00680 if (norf<ierr) // no inner iteration loop was required *** 00681 00682 { 00683 totl += dl; 00684 00685 if (j < 0) // no inner loop was required 00686 modif++; 00687 00688 if (modif>1) // two successive steps did not require inner loop 00689 { 00690 // length of arc can be increased 00691 dl*=2.0; 00692 if (dl>dlmax) 00693 dl=dlmax; 00694 modif=0; 00695 if (Mespr==1) 00696 fprintf(stdout,"\n arc-length is modified (dl increase) dl=%e",dl); 00697 } 00698 lambda+=dlambda; 00699 00700 // this part modifies parameter psi 00701 if (modif_psi) 00702 adapt_psi_coeff(i, ddlambda, ddlambda0, psi0, psi); 00703 00704 Mm->updateipval(); 00705 if (outres == yes) 00706 { 00707 compute_req_val(lcid); 00708 print_step(lcid, i+1, lambda, fa); 00709 print_flush(); 00710 } 00711 00712 // required value of lambda was attained for the load vector -> end of arclength 00713 if ((nlman->check_lambdar == on) && (check_rv_fp == off)) 00714 break; 00715 00716 // required total(cumulative) arc length has been attained => end of arclength 00717 if ((nlman->check_tot_al == on) && (check_max_tot_al==off)) 00718 break; 00719 } 00720 00721 00722 if (stop==1) // the equilibrium was attained 00723 { 00724 // this part modifies parameter psi 00725 if (modif_psi) 00726 adapt_psi_coeff(i, ddlambda, ddlambda0, psi0, psi); 00727 00728 if ((nlman->check_lambdar == on) && (check_rv_fp == off) && (fabs(lambda-nlman->lambdar) > nlman->errl)) 00729 check_rv_fp = on; 00730 } 00731 00732 if ((stop==0) || (stopj)) // stop == 0 the equilibrium was not attained in the prescribed number of iterations 00733 // stopj != 0 the dlambda could not be determined or divergence detected 00734 { 00735 if (dl == dlmin) // dl cannot be decreased 00736 break; 00737 // modification of the arc length 00738 dl/=2.0; 00739 if (dl<dlmin) 00740 dl=dlmin; 00741 00742 if (Mespr==1) 00743 { 00744 if (dl == dlmin) 00745 { 00746 fprintf (stdout,"\n arc length was decreased dl=dlmin=%e",dl); 00747 // fprintf (Out,"\n arc length was decreased dl=dlmin=%e",dl); 00748 } 00749 else 00750 { 00751 fprintf (stdout,"\n arc length was decreased dl=%e",dl); 00752 // fprintf (Out,"\n arc length was decreased dl=%e",dl); 00753 } 00754 } 00755 // restoring of left hand side vector 00756 copyv (r, ra, n); 00757 // restoring of lambda parameter 00758 lambda=blambda; 00759 } 00760 00761 // ------------------------------------ 00762 // finish of main iteration loop ---- 00763 // ------------------------------------ 00764 } 00765 00766 if ((i < ni-1) && (norf >= ierr)) 00767 { 00768 if (Mespr==1) 00769 { 00770 print_err("the step length was decreased to the minimum required value\n" 00771 " but the equilibrium could not be attained.\n" 00772 " FORCED output of the attained results was performed in this step.", 00773 __FILE__, __LINE__, __func__); 00774 } 00775 compute_req_val (lcid); 00776 print_step_forced(lcid, i+1, lambda, f); 00777 print_flush(); 00778 } 00779 00780 if (Mp->hdbcont.save_stat()) 00781 { 00782 // backup of the attained state is required 00783 if (Mespr==1) 00784 fprintf (stdout,"\n Creating backup file\n"); 00785 solver_save (ra, fp, i, lambda, dl, NULL, n); 00786 } 00787 00788 00789 delete [] r; 00790 delete [] fi; 00791 delete [] f; 00792 delete [] v; 00793 delete [] u; 00794 delete [] ddr; 00795 00796 return lambda; 00797 } 00798 */ 00799 00800 00801 /** 00802 Function performs calculation of one load/time step of the Newton-Raphson 00803 method for the given load case. Solved equation system does not contain 00804 time variable. 00805 00806 @param lcid - load case id 00807 @param nlman - pointer to structure conatining setup of the solver 00808 @param fa - attained load %vector 00809 @param ra - %vector of attained displacements 00810 @param fb - residual %vector or load %vector increment - righthand side %vector 00811 @param dr - %vector of displacement increment - lefthand side %vector 00812 @param fb - %vector of internal forces 00813 @param f - auxiliary %vector of actual load (call of solve_system rewrite it during solution) 00814 @param istep - time/load step id 00815 @param j - inner loop step id (output parameter, set to -1 if no inner loop was performed) 00816 @param li - initial value of time/load step id 00817 00818 @return The function returns reached lambda parameter. 00819 00820 Created by Tomas Koudelka, 11.2012 00821 */ 00822 /* 00823 double garclength_one_step(long lcid, nonlinman *nlman, double *fa, double *ra, double *fb, double *dr, double *fi, 00824 long istep, long &j, long li) 00825 { 00826 double *r; // displacements from previous step (backup copy) 00827 double *ddr; // vector of displacement increments 00828 double *u; // vector of corrections for computed displacements v 00829 double *fi; // vector of internal forces 00830 double *ra; // %vector of attained displacements 00831 double *fa; // attained load %vector 00832 double *fc; // constant part of load %vector (it is not influenced by lambda) 00833 double *fp; // proportional part of load %vector (it is increased by load coefficient lambda up to the nlman->lambdar) 00834 00835 double norfp; // norm of proportional load vector 00836 double psi; // proportional coefficient between load vector and displacement vector 00837 long modif_psi; // flag for automatic modification psi in dependence on the ddlambda 00838 double lambda; // actual value of load coefficient 00839 flagsw &check_rv_fp; // flag for checking of attainment of required value for load coefficient lambdar 00840 flagsw &check_max_tot_al; // flag fro checking maximu value of the length of arc 00841 00842 double *v; // vector of computed displacements for actual arclength step 00843 double *f; // auxiliary vector of actual load (call of solve_system rewrite it during solution) 00844 00845 double norv; // norm of displacement increments 00846 double norf; // norm of residual vector 00847 double norfa; // norm of attained load vector 00848 double dlambda; // load coefficient increment from the inner loop 00849 double ddlambda; // cummulative value of load coefficent increments ddlambda = \sum(dlambda) 00850 double ddlambda0;// initial/reference value of ddlambda 00851 double tmp; // auxiliary variable 00852 00853 // number of rows of the matrix 00854 n = Ndofm; 00855 // maximum number of iterations in one increment 00856 ini = nlman->niilal; 00857 // required error in inner loop 00858 ierr = nlman->erral; 00859 // maximum number of iterations in one increment 00860 ini = nlman->niilal; 00861 // required error in inner loop 00862 ierr = nlman->erral; 00863 // length of arc 00864 dl = nlman->dlal; 00865 // maximum length of arc 00866 dlmax = nlman->dlmaxal; 00867 // minimum length of arc 00868 dlmin = nlman->dlminal; 00869 00870 // Delta lambda - cumulative increment in one load increment 00871 // at the beginning of computation is always zero 00872 dlambda=0.0; 00873 00874 // assembling of the stiffness matrix 00875 sm_change = assemble_stiffness_matrix (lcid, istep, -1, li); 00876 00877 // solution of K(r).v=F 00878 Mp->ssle->solve_system(Gtm, Smat, v, f, Out); 00879 00880 // square of generalized norm of displacement increments 00881 norv = displincr(v,v,n); 00882 00883 if ((istep == 0) && (psi < 0.0) && (norfp > Mp->zero)) 00884 psi0 = psi = sqrt(norv/norfp/norfp); 00885 00886 // compute new dlambda increment 00887 tmp = norv+psi*psi*norfp*norfp; 00888 // compute new dlambda increment 00889 dlambda = dl/sqrt(tmp); 00890 00891 // check required value of load coefficient lambdar 00892 if ((check_rv_fp == on) && (dlambda+lambda > nlman->lambdar)) 00893 { 00894 check_rv_fp = off; 00895 dlambda = nlman->lambdar - lambda; 00896 dl = dlambda*sqrt(tmp); 00897 } 00898 // checking of maximum value of the total (cumulative) arc length 00899 if ((check_max_tot_al == on) && (dl+totl > nlman->max_tot_al)) 00900 { 00901 check_max_tot_al = off; 00902 dl = nlman->max_tot_al - totl; 00903 } 00904 00905 if (Mespr == 1) 00906 fprintf(Out, "Step %ld: ddlambda = %e, norv = %e, norfp = %e\n", i, ddlambda, norv, norfp); 00907 00908 // ddr[j]=dlambda*v[j]; 00909 copymultv(v, ddr, dlambda, n); 00910 00911 // ra[j]+=ddr[j]; 00912 addv(ra, ddr, n); 00913 00914 // update attained load vector 00915 // fa[j] = fc + (lambda+dlambda)*fp[j] is calculated 00916 addmultv(fc, fp, lambda+dlambda, fa, n); 00917 00918 00919 // initialize value of ddlambda before inner iteration loop 00920 ddlambda=dlambda; 00921 00922 // computation of internal forces 00923 internal_forces (lcid, fi); 00924 00925 // f[j] = fa[j]-fi[j] 00926 subv(fa, fi, f, n); 00927 00928 // computation of norms 00929 norf = normv(f, n); 00930 norfa = normv(fa, n); 00931 00932 // proportional norm is used in case of nonzero load vector 00933 // otherwise absolute norm is used for vector of unbalanced forces 00934 if (fabs(norfa) > Mp->zero) 00935 norf /= norfa; 00936 00937 if (norf<ierr) 00938 { 00939 // ****************************************** 00940 // no inner iteration loop is required *** 00941 // ****************************************** 00942 if (Mespr==1) 00943 fprintf(stdout,"\n increment %ld norf=%e ierr=%e ddlambda=%e dlambda=%e",i,norf,ierr, ddlambda, dlambda); 00944 00945 j = -1; 00946 return norf; 00947 } 00948 00949 // *************************************** 00950 // * inner iteration loop * 00951 // *************************************** 00952 if (Mespr==1) 00953 fprintf (stdout,"\n inner loop: 1 norf=%e ierr=%e ddlambda=%e dlambda=%e",norf,ierr,ddlambda, dlambda); 00954 stop =0; 00955 stopj=0; 00956 fillm(0.0, lsm_a); 00957 fillv(0.0, lsm_r); 00958 for (j=0; j<ini; j++) 00959 { //loop j 00960 // copy of actual step to problem description 00961 // it is used in connection with printing of output values 00962 Mp->jstep = j; 00963 00964 // eventual update of stiffness matrix 00965 Mespr = 0; // switch off the message printing for assembling of stiffness %matrix 00966 sm_change = assemble_stiffness_matrix (lcid,i,j,li); 00967 Mespr = 1; // switch on the message printing for assembling of stiffness %matrix 00968 00969 // solve displacement increments due to vector of unbalanced forces 00970 Mp->ssle->solve_system(Gtm, Smat, u, f, Out); 00971 if (sm_change) 00972 { 00973 // backup of the fp because the right hand side will be destroyed during the solution of the system of equations 00974 copyv(fp, f, n); 00975 Mp->ssle->solve_system(Gtm, Smat, v, f, Out); 00976 } 00977 00978 // backup vector of displacement increments from the previous step 00979 // f[j] = ddr[j] 00980 copyv(ddr, f, n); 00981 00982 // add computed correction of displacement increments due to unbalanced forces 00983 // ddr[j] = ddr[j]+u[j] 00984 addv(ddr, u, n); 00985 00986 // determination of increment of load parameter 00987 determine_dlambda (ddr, u, v, f, n, ddlambda,psi,norfp,dl,dlambda,stopj, nlman); 00988 if (stopj) 00989 break; 00990 00991 // check required value of load coefficient lambdar 00992 if ((nlman->check_lambdar == on) && (check_rv_fp == off) && 00993 (dlambda+ddlambda+lambda > nlman->lambdar)) 00994 dlambda = nlman->lambdar - ddlambda - lambda; 00995 00996 // actualize vector of displacements increments 00997 // ddr[j] += dlambda*v[j] 00998 addmultv(ddr, v, dlambda, n); 00999 01000 // actualize vector of attained displacements ra 01001 // ra[j] += u[j] + dlambda*v[j] 01002 addv(ra, u, n); 01003 addmultv(ra, v, dlambda, n); 01004 // update attained load vector 01005 // fa[k]+=dlambda*fp[k] is calculated 01006 addmultv(fa, fp, dlambda, n); 01007 //update_attained_load_vector(fa, fp, n, lambda+ddlambda, dlambda, nlman); 01008 01009 ddlambda+=dlambda; 01010 01011 01012 // computation of internal forces 01013 internal_forces (lcid,fi); 01014 // computation of residual forces f[k] = fa[k] - fi[k] 01015 subv(fa, fi, f, n); 01016 // norm of the vector of residuals 01017 norf=normv(f,n); 01018 // norm of the vector of attained forces 01019 norfa = normv(fa, n); 01020 if (norfa > zero) 01021 { 01022 // if the attained forces are nonzero, the norm of residuals is scaled by the norm of forces 01023 // this scaling is reasonable because nondimensional quantity is obtained 01024 norf /= norfa; 01025 } 01026 01027 if (Mespr==1) 01028 { 01029 fprintf (stdout,"\n increment %ld inner loop %ld norf=%e",i,j,norf); 01030 // fprintf (Out,"\n increment %ld inner loop %ld norf=%e",i,j,norf); 01031 } 01032 01033 // equilibrium was attained 01034 if (norf<ierr) 01035 { 01036 lambda+=ddlambda; 01037 Mm->updateipval(); 01038 compute_req_val (lcid); 01039 print_step(lcid, i+1, lambda, fa); 01040 print_flush(); 01041 totl += dl; 01042 stop=1; 01043 break; 01044 } 01045 01046 // divergence detection with help of least square method 01047 stopj = check_divergency(nlman, lsm_a, lsm_r, lsm_l, j, norf); 01048 if (stopj) 01049 break; 01050 01051 } // end of j 01052 } 01053 */ 01054 01055 01056 /** 01057 The function initializes dof numbers for the particular 01058 types of displacement control. It is used in the arclength method. 01059 01060 @return The function does not return anything, it performs initialization of 01061 Mp->nlman. 01062 01063 Created by JK, 01064 */ 01065 void seldofinit () 01066 { 01067 long i,j,k,l,ndofn; 01068 double x1,x2,y1,y2,z1,z2,length; 01069 01070 switch (Mp->nlman->displnorm){ 01071 case alldofs:{ break; } 01072 case seldofs: 01073 case seldofscoord:{ 01074 for (i=0;i<Mp->nlman->nsdofal;i++){ 01075 j=Mp->nlman->seldofal[i]; 01076 Mp->nlman->seldofal[i]=Mt->give_dof (Mp->nlman->selnodal[i],j)-1; 01077 } 01078 break; 01079 } 01080 case selecnodes:{ 01081 Mp->nlman->nsdofal=0; 01082 for (i=0;i<Mp->nlman->nsnal;i++){ 01083 Mp->nlman->nsdofal+=Mt->give_ndofn (Mp->nlman->selnodal[i]); 01084 } 01085 Mp->nlman->seldofal = new long [Mp->nlman->nsdofal]; 01086 k=0; 01087 for (i=0;i<Mp->nlman->nsnal;i++){ 01088 ndofn=Mt->give_ndofn (Mp->nlman->selnodal[i]); 01089 for (j=0;j<ndofn;j++){ 01090 l=Mt->give_dof (Mp->nlman->selnodal[i],j); 01091 if (l>0){ 01092 Mp->nlman->seldofal[k]=l-1; 01093 k++; 01094 } 01095 } 01096 } 01097 Mp->nlman->nsdofal=k; 01098 break; 01099 } 01100 case nodesdistincr:{ 01101 Mp->nlman->nsdofal=0; 01102 for (i=0;i<Mp->nlman->nsnal;i++){ 01103 Mp->nlman->nsdofal+=Mt->give_ndofn (Mp->nlman->selnodal[i]); 01104 } 01105 Mp->nlman->seldofal = new long [Mp->nlman->nsdofal]; 01106 k=0; 01107 for (i=0;i<Mp->nlman->nsnal;i++){ 01108 ndofn=Mt->give_ndofn (Mp->nlman->selnodal[i]); 01109 for (j=0;j<ndofn;j++){ 01110 Mp->nlman->seldofal[k]=Mt->give_dof (Mp->nlman->selnodal[i],j)-1; 01111 k++; 01112 } 01113 } 01114 x1=Gtm->gnodes[Mp->nlman->selnodal[0]].x; x2=Gtm->gnodes[Mp->nlman->selnodal[1]].x; 01115 y1=Gtm->gnodes[Mp->nlman->selnodal[0]].y; y2=Gtm->gnodes[Mp->nlman->selnodal[1]].y; 01116 z1=Gtm->gnodes[Mp->nlman->selnodal[0]].z; z2=Gtm->gnodes[Mp->nlman->selnodal[1]].z; 01117 length=sqrt((x2-x1)*(x2-x1)+(y2-y1)*(y2-y1)+(z2-z1)*(z2-z1)); 01118 Mp->nlman->nxal=(x2-x1)/length; Mp->nlman->nyal=(y2-y1)/length; Mp->nlman->nzal=(z2-z1)/length; 01119 break; 01120 } 01121 default:{ 01122 fprintf (stderr,"\n\n unknown displacement norm is required in function seldofinit (%s, line %d)",__FILE__,__LINE__); 01123 } 01124 } 01125 } 01126 01127 01128 01129 /** 01130 The function computes generalized norm of load %vector. 01131 It is used in the arclength method. 01132 01133 @param fc - array with components of constant load %vector 01134 @param n - the number of components of displincr array 01135 @param lambda - attained value of the load coefficient 01136 @param nlman - structure with setup of handling of lambda and lambda tresholds 01137 01138 @return The function returns required norm. 01139 01140 Created by JK, 01141 */ 01142 double loadincr (double *fp, long n, double lambda, nonlinman *nlman) 01143 { 01144 long i; 01145 double norm; 01146 double *loadincr = new double[n]; 01147 01148 norm = 0.0; 01149 01150 nullv(loadincr, n); 01151 if ((nlman->check_lambdar == off) || (lambda > nlman->lambdar)) 01152 addv(loadincr, fp, n); 01153 01154 switch (Mp->nlman->displnorm) 01155 { 01156 case alldofs: 01157 case nodesdistincr: 01158 norm = ss (loadincr,loadincr,n); 01159 break; 01160 case seldofs: 01161 case seldofscoord: 01162 case selecnodes: 01163 for (i=0;i<Mp->nlman->nsdofal;i++) 01164 norm+=loadincr[Mp->nlman->seldofal[i]]*loadincr[Mp->nlman->seldofal[i]]; 01165 break; 01166 default: 01167 break; 01168 } 01169 norm = sqrt(norm); 01170 delete [] loadincr; 01171 01172 return norm; 01173 } 01174 01175 01176 01177 /** 01178 The function computes square of generalized norm of %vector of 01179 displacemnt increments. The function returns either square of %vector norm computed 01180 from the selected componets or it returns increment of distance between 01181 selected nodes. It is used in the arclength method. 01182 01183 @param dr1 - array with components of the first %vector of displacement increments 01184 @param dr2 - array with components of the second %vector of displacement increments 01185 @param n - number of components of displincr array 01186 @param lcid - load case id 01187 01188 @return The function returns required norm. 01189 01190 Created by JK, 01191 */ 01192 double displincr (double *dr1, double *dr2, long n, long lcid) 01193 { 01194 long i; 01195 double norm,u1,u2,v1,v2,w1,w2,aux; 01196 01197 switch (Mp->nlman->displnorm){ 01198 case alldofs:{ 01199 norm = scprd(dr1,dr2,n); 01200 break; 01201 } 01202 case seldofs: 01203 case seldofscoord: 01204 case selecnodes:{ 01205 norm=0.0; 01206 for (i=0;i<Mp->nlman->nsdofal;i++){ 01207 if (Mp->nlman->seldofal[i]>0) 01208 norm+=dr1[Mp->nlman->seldofal[i]]*dr2[Mp->nlman->seldofal[i]]; 01209 else 01210 { 01211 aux=prdisplincr(i,lcid); 01212 norm += aux*aux; 01213 } 01214 } 01215 break; 01216 } 01217 case nodesdistincr:{ 01218 if (Mp->nlman->probdimal==2){ 01219 u1=0.0; u2=0.0; v1=0.0; v2=0.0; 01220 if (Mp->nlman->seldofal[0]>0) u1=dr1[Mp->nlman->seldofal[0]]; 01221 if (Mp->nlman->seldofal[1]>0) v1=dr1[Mp->nlman->seldofal[1]]; 01222 if (Mp->nlman->seldofal[2]>0) u2=dr2[Mp->nlman->seldofal[2]]; 01223 if (Mp->nlman->seldofal[3]>0) v2=dr2[Mp->nlman->seldofal[3]]; 01224 01225 norm=(u2-u1)*Mp->nlman->nxal + (v2-v1)*Mp->nlman->nyal; 01226 norm *= norm; 01227 } 01228 if (Mp->nlman->probdimal==3){ 01229 u1=0.0; u2=0.0; v1=0.0; v2=0.0; w1=0.0; w2=0.0; 01230 if (Mp->nlman->seldofal[0]>0) u1=dr1[Mp->nlman->seldofal[0]]; 01231 if (Mp->nlman->seldofal[1]>0) v1=dr1[Mp->nlman->seldofal[1]]; 01232 if (Mp->nlman->seldofal[2]>0) w1=dr1[Mp->nlman->seldofal[2]]; 01233 if (Mp->nlman->seldofal[3]>0) u2=dr2[Mp->nlman->seldofal[3]]; 01234 if (Mp->nlman->seldofal[4]>0) v2=dr2[Mp->nlman->seldofal[4]]; 01235 if (Mp->nlman->seldofal[5]>0) w2=dr2[Mp->nlman->seldofal[5]]; 01236 01237 norm=(u2-u1)*Mp->nlman->nxal + (v2-v1)*Mp->nlman->nyal + (w2-w1)*Mp->nlman->nzal; 01238 norm *= norm; 01239 } 01240 break; 01241 } 01242 default: 01243 print_err("unknown norm of displacement increment is required", __FILE__, __LINE__, __func__); 01244 } 01245 return norm; 01246 } 01247 01248 01249 01250 /** 01251 The function returns in the case of prescribed displacement in the i-th selected dof from nonlinman. 01252 01253 @param i - id of selected dof 01254 @param lcid - load case id 01255 01256 @return The function returns contribution of to the length of arc from prescribed displacement 01257 of i-th selected dof. 01258 01259 Created by Tomas Koudelka, 5.2014 01260 */ 01261 double prdisplincr(long i, long lcid) 01262 { 01263 double ret = 0.0; 01264 01265 if (Mp->nlman->seldofal[i] < 0) 01266 { 01267 if (Mb->lc[lcid].pd) 01268 // Mp->nlman->seldof has been already decreased by one => index of pd must be decreased by 2 01269 ret = Mb->lc[lcid].pd[0-Mp->nlman->seldofal[i]-2]; 01270 01271 if (Mb->lc[lcid+1].pd) 01272 // Mp->nlman->seldof has been already decreased by one => index of pd must be decreased by 2 01273 ret = Mb->lc[lcid+1].pd[0-Mp->nlman->seldofal[i]-2]; 01274 } 01275 return ret; 01276 } 01277 01278 01279 01280 /** 01281 Function computes coefficients of the quadratic equation used in the arclength method. 01282 The coefficients are computed with respect to type of load/displacement control. 01283 01284 @param ddr - increment of load coefficient lambda for the given step 01285 @param v - %vector of displacement increments 01286 @param n - total number of DOFs (i.e. size of ddr and v arrays) 01287 @param ddlambda - correction of increment of load coefficient lambda 01288 @param psi - proportional weight coeffient (swithces between load/displacement control) 01289 @param norfp - norm of proportional load %vector 01290 @param dl - length of arc 01291 @param a0 - zero order coefficient (output) 01292 @param a1 - first order coefficient (output) 01293 @param a2 - second order coefficient (output) 01294 @param lcid - load case id 01295 01296 @return The function returns evaluated coefficient in the parameters a0, a1 and a2. 01297 01298 Created by JK, 01299 */ 01300 void quadeqcoeff (double *ddr,double *v,long n,double ddlambda,double psi,double norfp,double dl, 01301 double &a0,double &a1,double &a2,long lcid) 01302 { 01303 double norddr, norv, norddrv; 01304 01305 switch (Mp->nlman->displnorm){ 01306 case alldofs: 01307 case seldofs: 01308 case seldofscoord: 01309 case selecnodes:{ 01310 a0 = norddr = displincr(ddr, ddr, n, lcid); 01311 a2 = displincr(v, v, n, lcid); 01312 norddrv = displincr(ddr, v, n, lcid); 01313 a1 = 2.0*norddrv; 01314 01315 a0 += ddlambda*ddlambda*psi*psi*norfp*norfp - dl*dl; 01316 a1 += 2.0*ddlambda*psi*psi*norfp*norfp; 01317 a2 += psi*psi*norfp*norfp; 01318 break; 01319 } 01320 case nodesdistincr:{ 01321 norddr = displincr(ddr, ddr, n, lcid); 01322 norv = displincr(v, v, n, lcid); 01323 norddrv = sqrt(norddr)*sqrt(norv); // it is possible to calculate it in this way because ddr and v are scalars in this case 01324 a0 = norddr + ddlambda*ddlambda*norfp*norfp*psi*psi-dl*dl; 01325 a1 = 2.0*(norddrv + ddlambda*norfp*norfp*psi*psi); 01326 a2 = norv + norfp*norfp*psi*psi; 01327 break; 01328 } 01329 default: 01330 print_err("unknown norm of displacement increment is required", __FILE__, __LINE__, __func__); 01331 } 01332 // quadeqcoef_log(); // debugging log 01333 } 01334 01335 01336 01337 /** 01338 The function determines increment of lambda (delta lambda). 01339 01340 @param ddr - Delta r + u 01341 @param u - u = K^{-1}*(fa-fi) 01342 @param v - v = K^{-1}*fp 01343 @param ddrprev - Delta r (not updated by the %vector u) 01344 @param n - number of unknowns 01345 @param ddlambda - Delta lambda 01346 @param psi - proportional coefficient between load and displacement vectors 01347 @param norfp - the norm of proportional %vector 01348 @param dl - length of the arc 01349 @param dlambda - delta lambda (output of this function) 01350 @param stop - indicator of end of iteration (output of this function) 01351 @param nlman - pointer to structure with arclength setup (it is used for full arclength method) 01352 @param lcid - load case id 01353 01354 @return The function returns delta lambda in the parameter dlambda and 01355 it sets the stop parameter to indicate the end of the arclength due 01356 to nonsolvable constrained conditions. 01357 01358 Created by JK, 30.8.2010 01359 Modified 01360 */ 01361 void determine_dlambda (double *ddr, double *u, double *v, double *ddrprev, long n, 01362 double ddlambda, double psi, double norfp, double dl, 01363 double &dlambda, long &stop, nonlinman *nlman, long lcid) 01364 { 01365 long k,numr; 01366 double l1,l2,aux,a0,a1,a2,ss1,ss2,ss3,ss4,ss5,nom,denom; 01367 01368 switch (Mp->nlman->dlam) 01369 { 01370 case nodetermination:{ 01371 print_err("some reasonable type of lambda determination has to be selected",__FILE__,__LINE__,__func__); 01372 break; 01373 } 01374 case minvalue: 01375 case maxvalue: 01376 case minangle:{ 01377 // coefficient of quadratic equation 01378 quadeqcoeff (ddr,v,n,ddlambda,psi,norfp,dl,a0,a1,a2,lcid); 01379 // solution of quadratic equation 01380 numr = solv_polynom_2(a2, a1, a0, l1, l2); 01381 switch (numr){ 01382 case -1:{ 01383 print_err("infinite number of solution of constrained condition\n" 01384 "(all coefficients of quadratic equation are equal to zero", __FILE__, __LINE__, __func__); 01385 abort(); 01386 break; 01387 } 01388 case 0:{ 01389 print_err("nonsolvable constrained condition in function arclength", __FILE__, __LINE__, __func__); 01390 stop=0; 01391 break; 01392 } 01393 case 1:{ 01394 dlambda = l1; 01395 break; 01396 } 01397 default:{} 01398 } 01399 break; 01400 } 01401 case linearizedmeth:{ 01402 break; 01403 } 01404 case fullmethod: 01405 // unbalanced length of arc 01406 // a0 = ddrprev*ddrprev + ddlambda^2 * psi^2 * fp * fp - dl^2 01407 a0 = displincr(ddrprev, ddrprev, n, lcid); 01408 a0 += ddlambda*ddlambda*psi*psi*norfp*norfp - dl*dl; 01409 break; 01410 default:{ 01411 print_err("unknown type of detemination of lambda parameter is required",__FILE__,__LINE__,__func__); 01412 } 01413 } 01414 01415 01416 switch (Mp->nlman->dlam){ 01417 case nodetermination:{ 01418 print_err("some reasonable type of lambda determination has to be selected",__FILE__,__LINE__,__func__); 01419 break; 01420 } 01421 case minvalue:{ 01422 if (fabs(l1)>fabs(l2)) 01423 dlambda=l2; 01424 else 01425 dlambda=l1; 01426 break; 01427 } 01428 case maxvalue:{ 01429 if (fabs(l1)>fabs(l2)) 01430 dlambda=l1; 01431 else 01432 dlambda=l2; 01433 break; 01434 } 01435 case minangle:{ 01436 ss1=0.0; ss2=0.0; 01437 ss3=0.0; ss4=0.0; ss5=0.0; 01438 for (k=0;k<n;k++){ 01439 ss1+=(ddr[k]+l1*v[k])*ddrprev[k]; 01440 ss2+=(ddr[k]+l2*v[k])*ddrprev[k]; 01441 ss3+=(ddr[k]+l1*v[k])*(ddr[k]+l1*v[k]); 01442 ss4+=(ddr[k]+l2*v[k])*(ddr[k]+l2*v[k]); 01443 ss5+=ddrprev[k]*ddrprev[k]; 01444 } 01445 if (ss1/sqrt(ss3)/sqrt(ss5)>ss2/sqrt(ss4)/sqrt(ss5)) 01446 dlambda=l1; 01447 else 01448 dlambda=l2; 01449 break; 01450 } 01451 case linearizedmeth:{ 01452 nom=0.0; 01453 denom=0.0; 01454 for (k=0;k<n;k++){ 01455 nom-=ddrprev[k]*(ddr[k]-ddrprev[k]); 01456 denom+=ddrprev[k]*v[k]; 01457 } 01458 denom+=psi*psi*norfp*norfp*ddlambda; 01459 01460 dlambda=nom/denom; 01461 break; 01462 } 01463 case fullmethod: 01464 // - a0 - ddrprev*K^{-1}*(fa-fi) 01465 // dlambda = ---------------------------------------------- 01466 // 2*(ddrprev*K^{-1}*fp + ddlambda*psi^2 * fp*fp) 01467 aux = displincr(ddrprev, u, n, lcid); 01468 dlambda = -a0 - aux; 01469 aux = displincr(ddrprev, v, n, lcid); 01470 dlambda /= 2.0*(aux + ddlambda*psi*psi*norfp*norfp); 01471 break; 01472 default:{ 01473 print_err("unknown type of detemination of lambda parameter is required",__FILE__,__LINE__,__func__); 01474 } 01475 } 01476 } 01477 01478 01479 01480 /** 01481 The function modifies the psi coefficient according to the 01482 ddlambda increment adaptively. If the increment decreases comparing to 01483 the initial/reference ddlambda increment (ddlamda0) than the value should be also decreased 01484 and thus the arclength is controlled rather by the displacements than by the load. 01485 The initial/reference value of the load increment ddlambda0 and proportional coefficient psi0 01486 are updated at the beginning of the main iteration loop and if the ddlambda changes 01487 the sign comparing to the ddlambda 0 (peak point of loading curve). 01488 01489 @param i - step id 01490 @param ddlambda - actual increment of the load coefficeint lambda 01491 @param ddlambda0 - initial/reference value of the load coefficeint lambda increment (input/output parameter) 01492 @param psi0 - initial/reference value of the proportional coefficeint (input/output parameter) 01493 @param psi - actual value of the proportional coefficeint (input/output parameter) 01494 01495 @return The function returns the actual value of psi in the correspondign 01496 parameter. Also the parameters ddlambda0 and psi0 can be changed. 01497 01498 Created by Tomas Koudelka 07.2011 01499 */ 01500 void adapt_psi_coeff(long i, double ddlambda, double &ddlambda0, double &psi0, double &psi) 01501 { 01502 // at this point, ddlambda = dlambda 01503 // store initial values of load coefficient increment 01504 if (i==0) 01505 ddlambda0 = ddlambda; 01506 01507 // load coefficient increment changed sign -> peak was attained -> 01508 // actualize values of stored initail load coefficient increment and initial proportional coefficient 01509 if ((sgn(ddlambda0) != sgn(ddlambda)) && (sgn(ddlambda) != 0.0)) 01510 { 01511 psi0= psi; 01512 ddlambda0 = ddlambda; 01513 } 01514 01515 // automatic correction of proportional coefficient 01516 // with respect to actual and initial load coefficient increment 01517 if (ddlambda0 != 0.0) 01518 psi = psi0*fabs(ddlambda/ddlambda0); 01519 01520 return; 01521 }