00001 #include "nfssolver.h" 00002 #include "nssolver.h" 00003 #include "global.h" 00004 #include "globmat.h" 00005 #include "loadcase.h" 00006 #include "gmatrix.h" 00007 #include "mechprint.h" 00008 #include "flsubdom.h" 00009 #include "mathem.h" 00010 00011 00012 /** 00013 function solves problem with floating subdomains 00014 incremental analysis is used, compliance can vary 00015 00016 it may be used for pull-out tests 00017 00018 JK, 23. 11. 2012 00019 */ 00020 void solve_incremental_floating_subdomains () 00021 { 00022 long i,n,nm,lcid,nincr; 00023 double proppar; 00024 double *lhs,*rhs,*d,*lambda; 00025 00026 // load case id, must be zero 00027 lcid=0; 00028 // the number of increments 00029 nincr=Mp->nincr; 00030 // the number of primal DOFs 00031 n=Ndofm; 00032 // the number of dual DOFs - Lagrange multipliers 00033 nm=Mp->ssle->feti.ndofcp; 00034 // proportional parameter 00035 proppar=1.0/nincr; 00036 00037 // array for total displacements 00038 d = new double [n]; 00039 fillv (0.0,d,n); 00040 // array for total Lagrange multipliers 00041 lambda = new double [nm]; 00042 fillv (0.0,lambda,nm); 00043 // array for increments of displacements, array is allocated in Lsrs 00044 lhs = Lsrs->give_lhs (lcid); 00045 // array for increments of load, array is allocated in Lsrs 00046 rhs = Lsrs->give_rhs (lcid); 00047 00048 // assembling of the right hand side 00049 mefel_right_hand_side (lcid,rhs); 00050 // computation of load increment 00051 cmulv (proppar,rhs,n); 00052 00053 // assembling of the stiffness matrix 00054 stiffness_matrix (lcid); 00055 00056 Mp->ssle->feti.matrices_assembl (Smat,Out); 00057 Mp->ssle->feti.vectors_assembl (rhs,Out); 00058 00059 print_init(-1, "wt"); 00060 00061 // loop over the number of increments 00062 for (i=0;i<nincr;i++){ 00063 00064 Mp->ssle->feti.define_b (Out); 00065 Mp->ssle->feti.define_h (Out); 00066 00067 // solution of equation system 00068 Mp->ssle->feti.mpcg (Out); 00069 Mp->ssle->feti.nodalunknowns (d,Out); 00070 00071 // total displacements are computed 00072 addv (lhs,d,n); 00073 00074 print_init(-1, "wt"); 00075 // computes and prints required quantities 00076 compute_req_val (lcid); 00077 print_step(lcid, i+1, 0.0, NULL); 00078 00079 }// end of the loop over the number of increments 00080 00081 for (i=0;i<Lsrs->nlc;i++){ 00082 // computes and prints required quantities 00083 //compute_ipstrains (i); 00084 //compute_ipstresses (i); 00085 compute_req_val (i); 00086 print_step(i, 0, 0.0, NULL); 00087 } 00088 print_close(); 00089 00090 } 00091 00092 00093 /** 00094 function solves linear problem with floating subdomains 00095 it may be used for pull-out tests 00096 00097 JK, 5.11.2005 00098 */ 00099 /* 00100 00101 void solve_incremental_floating_subdomains () 00102 { 00103 long i,lcid,n,nincr; 00104 double *d,*lhs,*rhs; 00105 00106 // load case id, must be zero 00107 lcid=0; 00108 // the number of unknowns 00109 n=Ndofm; 00110 // the number of increments 00111 nincr=Mp->nincr; 00112 00113 Fsd = new flsubdom; 00114 00115 // array of nodal displacements 00116 lhs=Lsrs->give_lhs (lcid); 00117 // array of nodal forcess 00118 rhs=Lsrs->give_rhs (lcid); 00119 // array of increments of nodal displacements 00120 d = new double [n]; 00121 00122 // assembling of the right hand side 00123 mefel_right_hand_side (lcid,rhs); 00124 00125 // assembling of the stiffness matrix 00126 stiffness_matrix (lcid); 00127 00128 Fsd->initialization (Ndofm,Mp->ense,rhs); 00129 00130 print_init(-1, "wt"); 00131 for (i=0;i<nincr;i++){ 00132 00133 Fsd->solve_lin_alg_system (lhs,rhs); 00134 00135 addv (lhs,d,n); 00136 00137 Fsd->add_mult (1.0); 00138 Fsd->mult_correction (); 00139 00140 for (i=0;i<Lsrs->nlc;i++){ 00141 //compute_ipstrains (i); 00142 //compute_ipstresses (i); 00143 compute_req_val (i); 00144 print_step(i, 0, 0.0, NULL); 00145 } 00146 00147 } 00148 00149 print_close(); 00150 00151 00152 } 00153 */ 00154 00155 /** 00156 function solves nonlinear problem of floating subdomain 00157 the nonlinearity is hidden in contact stresses and in 00158 material models, like plasticity, etc. 00159 00160 function solves the problem by the arc-length method (continuation method) 00161 this function is modification of the function arclength 00162 from the file nssolver.cpp 00163 00164 d stands for delta 00165 dd stands for capital delta 00166 00167 fc - nonproportional %vector 00168 fp - proportional %vector 00169 n - number of unknowns 00170 00171 JK, 24.6.2006 00172 */ 00173 /* 00174 void solve_nonlinear_floating_subdomains () 00175 { 00176 long i,j,k,n,l,ni,ini,stop,modif,li,newmesh,numr,lcid,nlm; 00177 double a0,a1,a2,l1,l2,dl,dlmax,dlmin,psi,lambda,blambda,dlambda,ddlambda,norf,norfp,norv,norfa,zero,ierr; 00178 double lambdao,ss1,ss2,ss3,ss4,ss5; 00179 double *r,*ra,*ddr,*u,*v,*f,*fa,*fp,*fc,*fi,*rhs,*blm; 00180 00181 // allocation of object floating subdomains 00182 Fsd = new flsubdom; 00183 00184 // only one load case can be solved 00185 // it is nonlinear computation and the superposition is useless 00186 lcid=0; 00187 00188 // number of rows of the matrix 00189 n = Ndofm; 00190 // maximum number of increments 00191 ni = Mp->nlman->nial; 00192 // maximum number of iterations in one increment 00193 ini = Mp->nlman->niilal; 00194 // computer zero 00195 zero = Mp->zero; 00196 // required error in inner loop 00197 ierr = Mp->nlman->erral; 00198 // length of arc 00199 dl = Mp->nlman->dlal; 00200 // maximum length of arc 00201 dlmax = Mp->nlman->dlmaxal; 00202 // minimum length of arc 00203 dlmin = Mp->nlman->dlminal; 00204 // displacement-loading driving switch 00205 psi = Mp->nlman->psial; 00206 00207 // array for right hand side 00208 rhs=Lsrs->give_rhs (0); 00209 // assembling of right hand side (load vector) 00210 mefel_right_hand_side (lcid,rhs); 00211 00212 dlambda=0.0; 00213 00214 // allocation of auxiliary arrays 00215 r = new double [n]; 00216 ddr = new double [n]; 00217 u = new double [n]; 00218 v = new double [n]; 00219 f = new double [n]; 00220 fa = new double [n]; 00221 fi = new double [n]; 00222 00223 // initialization phase 00224 // vector of nodal displacements 00225 ra = Lsrs->give_lhs (lcid); 00226 // vector of proportional load 00227 fp = Lsrs->give_rhs (lcid*2); 00228 // vector of nonproportional load 00229 fc = Lsrs->give_rhs (lcid*2+1); 00230 00231 00232 if (Mp->nlman->hdbackupal == 2){ // termit 00233 arclopen (li,n,lambda,dl,ra,fp); 00234 } 00235 else{ 00236 lambda=0.0; 00237 lambdao=0.0; 00238 li=0; 00239 } 00240 00241 00242 // norm of proportionality vector 00243 norfp = ss(fp,fp,n); 00244 modif=0; 00245 00246 00247 // assembling reached load vector 00248 for (j=0;j<n;j++){ 00249 fa[j]=fc[j]+(lambda+dlambda)*fp[j]; 00250 } 00251 if (li) 00252 print_init(-1, "at"); 00253 else 00254 { 00255 print_init(-1, "wt"); 00256 print_step(lcid, li, lambda, fa); 00257 } 00258 00259 00260 // assembling of stiffness matrix 00261 stiffness_matrix (lcid); 00262 00263 // initialization of FETI method 00264 // determination of number of Lagrange multipliers 00265 // computation of rigid body motions 00266 Fsd->initialization (Ndofm,Mp->ense,f); 00267 00268 // number of Lagrange multipliers 00269 nlm = Fsd->nlm; 00270 00271 00272 00273 // *************************** 00274 // main iteration loop **** 00275 // *************************** 00276 for (i=li;i<ni;i++){ 00277 00278 fprintf (Out,"\n\n arc-length prirustek %ld",i); 00279 00280 stop=1; // termit 00281 // backup of left hand side vector 00282 copyv (ra, r, n); 00283 // backup of reached lambda parameter 00284 blambda=lambda; 00285 // backup of reached Lagrange multipliers 00286 copyv (Fsd->tw,blm,nlm); 00287 00288 //fprintf (stdout,"\n\n arc-length: increment %ld lambda %e dl %e",i,lambda,dl); 00289 // if (Mespr==1){ 00290 // fprintf (Out,"\n\n *******************************************************************"); 00291 // fprintf (Out,"\n arc-length: increment %ld, lambda=%e dl=%e",i,lambda,dl); 00292 // } 00293 00294 00295 // assembling of tangent stiffness matrix 00296 if ((Mp->nlman->stmat==tangent_stiff) || (i == li)){ 00297 //stiffness_matrix (lcid); 00298 //Fsd->initialization (Ndofm,Mp->ense,f); 00299 } 00300 00301 stiffness_matrix (lcid); 00302 00303 00304 // backup of the fp, in ldl_sky the right hand side will be destroyed 00305 copyv(fp, f, n); 00306 00307 Fsd->initialization (Ndofm,Mp->ense,f); 00308 00309 Fsd->solve_lin_alg_system (v,f); 00310 00311 // solution of K(r).v=F 00312 //Smat->solve_system (v,f); 00313 00314 // generalized norm of displacement increments 00315 norv = displincr (v,n); 00316 00317 // compute new dlambda increment 00318 dlambda = dl/sqrt(norv+psi*psi*norfp); 00319 00320 Fsd->add_mult (dlambda); 00321 Fsd->mult_correction (); 00322 00323 // novinka 00324 if (Fsd->nch!=0){ 00325 for (j=0;j<ini;j++){ 00326 00327 fprintf (stdout,"\n cyklime v casti jedna %ld",j); 00328 00329 // restoring of Lagrange multipliers 00330 copyv (blm,Fsd->tw,nlm); 00331 00332 Fsd->solve_lin_alg_system (v,f); 00333 00334 // generalized norm of displacement increments 00335 norv = displincr (v,n); 00336 00337 // compute new dlambda increment 00338 dlambda = dl/sqrt(norv+psi*psi*norfp); 00339 00340 Fsd->add_mult (dlambda); 00341 Fsd->mult_correction (); 00342 00343 if (Fsd->nch==0) 00344 break; 00345 } 00346 } 00347 // konec novinky 00348 00349 00350 for (j=0;j<n;j++){ 00351 ddr[j]=dlambda*v[j]; 00352 ra[j]+=ddr[j]; 00353 fa[j]=fc[j]+(lambda+dlambda)*fp[j]; 00354 } 00355 ddlambda=dlambda; 00356 00357 // computation of internal forces 00358 internal_forces (lcid,fi); 00359 subv(fa, fi, f, n); 00360 norf=sizev(f,n); 00361 norfa = sizev(fa, n); 00362 //if (norfa<1.0e-8){} 00363 //else norf /= norfa; 00364 norf /= norfa; 00365 00366 00367 //if (Mespr==1) fprintf (stdout,"\n ddlambda %e dlambda %e norf %e ierr %e",ddlambda,dlambda,norf,ierr); 00368 if (Mespr==1) fprintf (stdout,"\n increment %ld norv %e norf %e",i,norv,norf); 00369 00370 if (norf<ierr){ 00371 // ****************************************** 00372 // no inner iteration loop is required *** 00373 // ****************************************** 00374 modif++; 00375 00376 if (modif>1){ 00377 // arc length modification 00378 dl*=2.0; 00379 if (dl>dlmax) 00380 dl=dlmax; 00381 modif=0; 00382 if (Mespr==1){ 00383 fprintf (stdout,"\n arc-length must be modified (dl increase) dl=%e",dl); 00384 // fprintf (Out,"\n arc-length must be modified (dl increase) dl=%e",dl); 00385 } 00386 } 00387 lambda+=dlambda; 00388 Mm->updateipval (); 00389 compute_req_val (lcid); 00390 print_step(lcid, i+1, lambda, fa); 00391 print_flush(); 00392 } 00393 else{ 00394 // **************************** 00395 // inner iteration loop **** 00396 // **************************** 00397 stop=0; 00398 for (j=0;j<ini;j++){ 00399 if (Mp->nlman->stmat==tangent_stiff){ 00400 stiffness_matrix (lcid); 00401 Fsd->initialization (Ndofm,Mp->ense,f); 00402 } 00403 00404 fprintf (Out,"\n arc-length vnitrni smycka %ld %ld",i,j); 00405 00406 00407 copyv(f, buf, n); 00408 00409 Fsd->solve_lin_alg_system (u,f); 00410 00411 // back substitution 00412 //Smat->solve_system (u,f); 00413 00414 // backup of the vector ddr 00415 copyv(ddr, buddr, n); 00416 00417 00418 copyv(ddr, f, n); 00419 addv(ddr, u, n); 00420 // coefficient of quadratic equation 00421 quadeqcoeff (ddr,v,n,ddlambda,psi,norfp,dl,a0,a1,a2); 00422 00423 // solution of quadratic equation 00424 numr = solv_polynom_2(a2, a1, a0, l1, l2); 00425 switch (numr) 00426 { 00427 case -1 : 00428 fprintf (stderr,"\n\n infinite number of solution of constrained condition in function arclength"); 00429 break; 00430 case 0 : 00431 fprintf (stderr,"\n\n nonsolvable constrained condition in function arclength"); 00432 break; 00433 case 1 : 00434 dlambda = l1; 00435 break; 00436 default : 00437 break; 00438 } 00439 ss1=0.0; ss2=0.0; 00440 ss3=0.0; ss4=0.0; ss5=0.0; 00441 for (k=0;k<n;k++){ 00442 ss1+=(ddr[k]+l1*v[k])*f[k]; 00443 ss2+=(ddr[k]+l2*v[k])*f[k]; 00444 ss3+=(ddr[k]+l1*v[k])*(ddr[k]+l1*v[k]); 00445 ss4+=(ddr[k]+l2*v[k])*(ddr[k]+l2*v[k]); 00446 ss5+=f[k]*f[k]; 00447 } 00448 00449 if (fabs(l1)>fabs(l2)) dlambda=l2; 00450 else dlambda=l1; 00451 00452 copyv (Fsd->tw,bulm,nlm); 00453 00454 Fsd->add_mult (dlambda); 00455 Fsd->mult_correction (); 00456 00457 // novinka 00458 if (Fsd->nch!=0){ 00459 for (l=0;l<ini;l++){ 00460 00461 fprintf (stdout,"\n cyklime v casti dve %ld",j); 00462 00463 00464 // restoring of Lagrange multipliers 00465 copyv (bulm,Fsd->tw,nlm); 00466 copyv (buddr,ddr,n); 00467 copyv (buf,f,n); 00468 00469 Fsd->solve_lin_alg_system (u,f); 00470 00471 copyv(ddr, f, n); 00472 addv(ddr, u, n); 00473 // coefficient of quadratic equation 00474 quadeqcoeff (ddr,v,n,ddlambda,psi,norfp,dl,a0,a1,a2); 00475 00476 // solution of quadratic equation 00477 numr = solv_polynom_2(a2, a1, a0, l1, l2); 00478 switch (numr) 00479 { 00480 case -1 : 00481 fprintf (stderr,"\n\n infinite number of solution of constrained condition in function arclength"); 00482 break; 00483 case 0 : 00484 fprintf (stderr,"\n\n nonsolvable constrained condition in function arclength"); 00485 break; 00486 case 1 : 00487 dlambda = l1; 00488 break; 00489 default : 00490 break; 00491 } 00492 ss1=0.0; ss2=0.0; 00493 ss3=0.0; ss4=0.0; ss5=0.0; 00494 for (k=0;k<n;k++){ 00495 ss1+=(ddr[k]+l1*v[k])*f[k]; 00496 ss2+=(ddr[k]+l2*v[k])*f[k]; 00497 ss3+=(ddr[k]+l1*v[k])*(ddr[k]+l1*v[k]); 00498 ss4+=(ddr[k]+l2*v[k])*(ddr[k]+l2*v[k]); 00499 ss5+=f[k]*f[k]; 00500 } 00501 00502 if (fabs(l1)>fabs(l2)) dlambda=l2; 00503 else dlambda=l1; 00504 00505 00506 Fsd->add_mult (dlambda); 00507 Fsd->mult_correction (); 00508 00509 if (Fsd->nch==0) 00510 break; 00511 } 00512 } 00513 // konec novinky 00514 00515 00516 00517 00518 for (k=0;k<n;k++){ 00519 ddr[k]+=dlambda*v[k]; 00520 ra[k]+=u[k]+dlambda*v[k]; 00521 fa[k]+=dlambda*fp[k]; 00522 } 00523 ddlambda+=dlambda; 00524 00525 //fprintf (stdout," ddlambda %e",ddlambda); 00526 00527 //fprintf (Out,"\n ddlambda %e dlambda %e",ddlambda,dlambda); 00528 // computation of internal forces 00529 internal_forces (lcid,fi); 00530 subv(fa, fi, f, n); 00531 norf=sizev(f,n); 00532 norfa = sizev(fa, n); 00533 norf /= norfa; 00534 //if (norfa<1.0e-8){} 00535 //else norf /= norfa; 00536 00537 if (Mespr==1){ 00538 fprintf (stdout,"\n norf=%e ierr=%e",norf,ierr); 00539 //fprintf (Out,"\n increment %ld inner loop %ld norf=%e",i,j,norf); 00540 } 00541 00542 if (norf<ierr){ 00543 lambda+=ddlambda; 00544 Mm->updateipval (); 00545 stop=1; 00546 compute_req_val (lcid); 00547 print_step(lcid, i+1, lambda, fa); 00548 print_flush(); 00549 break; 00550 } 00551 } 00552 modif=0; 00553 if (stop==0){ 00554 // modification of the arc length 00555 dl/=2.0; 00556 if (dl<dlmin){ 00557 dl=dlmin; 00558 break; 00559 } 00560 if (Mespr==1){ 00561 fprintf (stdout,"\n arc length must be modified (dl decrease) dl=%e",dl); 00562 //fprintf (Out,"\n arc length must be modified (dl decrease) dl=%e",dl); 00563 } 00564 // restoring of left hand side vector 00565 copyv(r, ra, n); 00566 // restoring of lambda parameter 00567 lambda=blambda; 00568 // restoring of Lagrange multipliers 00569 copyv (blm,Fsd->tw,nlm); 00570 } 00571 } 00572 00573 fprintf (stdout,"\n increment %ld total lambda %e",i,lambda); 00574 00575 if (stop==0) 00576 continue; 00577 00578 } 00579 00580 // ------------------------------------ 00581 // finish of main iteration loop ---- 00582 // ------------------------------------ 00583 00584 00585 fprintf (stdout,"\n\n multiplikatory %lf",Fsd->tw[0]); 00586 00587 if (Mp->nlman->hdbackupal==1) 00588 arclsave (i,n,lambda,dl,ra,fp); 00589 00590 print_close(); 00591 delete [] r; 00592 delete [] fi; 00593 delete [] fa; 00594 delete [] f; 00595 delete [] v; 00596 delete [] u; 00597 delete [] ddr; 00598 00599 } 00600 */ 00601 00602 00603 /** 00604 function solves nonlinear problem of floating subdomain 00605 the nonlinearity is hidden in contact stresses and in 00606 material models, like plasticity, etc. 00607 00608 function solves the problem by the arc-length method (continuation method) 00609 this function is modification of the function arclength 00610 from the file nssolver.cpp 00611 00612 d stands for delta 00613 dd stands for capital delta 00614 00615 fc - nonproportional %vector 00616 fp - proportional %vector 00617 n - number of unknowns 00618 00619 JK, 20.6.2006 (25.7.2001) 00620 */ 00621 /* 00622 void solve_nonlinear_floating_subdomains_24_6 () 00623 { 00624 long i,j,k,n,l,ni,ini,stop,modif,li,numr,lcid,nlm; 00625 double a0,a1,a2,l1,l2,dl,dlmax,dlmin,psi,lambda,blambda,dlambda,ddlambda,norf,norfp,norv,norfa,zero,ierr; 00626 double lambdao,ss1,ss2,ss3,ss4,ss5; 00627 double *r,*ra,*ddr,*u,*v,*f,*fa,*fp,*fc,*fi,*rhs,*blm; 00628 00629 Fsd = new flsubdom; 00630 00631 lcid=0; 00632 00633 rhs=Lsrs->give_rhs (0); 00634 mefel_right_hand_side (lcid,rhs); 00635 00636 // number of rows of the matrix 00637 n = Ndofm; 00638 // maximum number of increments 00639 ni = Mp->nlman->nial; 00640 // maximum number of iterations in one increment 00641 ini = Mp->nlman->niilal; 00642 // computer zero 00643 zero = Mp->zero; 00644 // required error in inner loop 00645 ierr = Mp->nlman->erral; 00646 // length of arc 00647 dl = Mp->nlman->dlal; 00648 // maximum length of arc 00649 dlmax = Mp->nlman->dlmaxal; 00650 // minimum length of arc 00651 dlmin = Mp->nlman->dlminal; 00652 // displacement-loading driving switch 00653 psi = Mp->nlman->psial; 00654 00655 dlambda=0.0; 00656 00657 // allocation of auxiliary arrays 00658 r = new double [n]; 00659 ddr = new double [n]; 00660 u = new double [n]; 00661 v = new double [n]; 00662 f = new double [n]; 00663 fa = new double [n]; 00664 fi = new double [n]; 00665 00666 // initialization phase 00667 // vector of nodal displacements 00668 ra = Lsrs->give_lhs (lcid); 00669 // vector of proportional load 00670 fp = Lsrs->give_rhs (lcid*2); 00671 // vector of nonproportional load 00672 fc = Lsrs->give_rhs (lcid*2+1); 00673 00674 00675 if (Mp->nlman->hdbackupal == 2){ // termit 00676 arclopen (li,n,lambda,dl,ra,fp); 00677 } 00678 else{ 00679 lambda=0.0; 00680 lambdao=0.0; 00681 li=0; 00682 } 00683 00684 00685 // norm of proportionality vector 00686 norfp = ss(fp,fp,n); 00687 modif=0; 00688 00689 00690 // assembling reached load vector 00691 for (j=0;j<n;j++){ 00692 fa[j]=fc[j]+(lambda+dlambda)*fp[j]; 00693 } 00694 if (li) 00695 print_init(-1, "at"); 00696 else 00697 { 00698 print_init(-1, "wt"); 00699 print_step(lcid, li, lambda, fa); 00700 } 00701 00702 00703 00704 stiffness_matrix (lcid); 00705 Fsd->initialization (Ndofm,Mp->ense,f); 00706 // number of Lagrange multipliers 00707 nlm = Fsd->nlm; 00708 00709 blm = new double [nlm]; 00710 00711 double *buddr; 00712 buddr = new double [n]; 00713 double *buf; 00714 buf = new double [n]; 00715 double *bulm; 00716 bulm = new double [nlm]; 00717 00718 // *************************** 00719 // main iteration loop **** 00720 // *************************** 00721 for (i=li;i<ni;i++){ 00722 00723 fprintf (Out,"\n\n arc-length prirustek %ld",i); 00724 00725 stop=1; // termit 00726 // backup of left hand side vector 00727 copyv (ra, r, n); 00728 // backup of reached lambda parameter 00729 blambda=lambda; 00730 // backup of reached Lagrange multipliers 00731 copyv (Fsd->tw,blm,nlm); 00732 00733 //fprintf (stdout,"\n\n arc-length: increment %ld lambda %e dl %e",i,lambda,dl); 00734 if (Mespr==1){ 00735 fprintf (Out,"\n\n *******************************************************************"); 00736 fprintf (Out,"\n arc-length: increment %ld, lambda=%e dl=%e",i,lambda,dl); 00737 } 00738 00739 00740 // assembling of tangent stiffness matrix 00741 if ((Mp->nlman->stmat==tangent_stiff) || (i == li)){ 00742 //stiffness_matrix (lcid); 00743 //Fsd->initialization (Ndofm,Mp->ense,f); 00744 } 00745 00746 stiffness_matrix (lcid); 00747 00748 00749 // backup of the fp, in ldl_sky the right hand side will be destroyed 00750 copyv(fp, f, n); 00751 00752 Fsd->initialization (Ndofm,Mp->ense,f); 00753 00754 Fsd->solve_lin_alg_system (v,f); 00755 00756 // solution of K(r).v=F 00757 //Smat->solve_system (v,f); 00758 00759 // generalized norm of displacement increments 00760 norv = displincr (v,n); 00761 00762 // compute new dlambda increment 00763 dlambda = dl/sqrt(norv+psi*psi*norfp); 00764 00765 Fsd->add_mult (dlambda); 00766 Fsd->mult_correction (); 00767 00768 // novinka 00769 if (Fsd->nch!=0){ 00770 for (j=0;j<ini;j++){ 00771 00772 fprintf (stdout,"\n cyklime v casti jedna %ld",j); 00773 00774 // restoring of Lagrange multipliers 00775 copyv (blm,Fsd->tw,nlm); 00776 00777 Fsd->solve_lin_alg_system (v,f); 00778 00779 // generalized norm of displacement increments 00780 norv = displincr (v,n); 00781 00782 // compute new dlambda increment 00783 dlambda = dl/sqrt(norv+psi*psi*norfp); 00784 00785 Fsd->add_mult (dlambda); 00786 Fsd->mult_correction (); 00787 00788 if (Fsd->nch==0) 00789 break; 00790 } 00791 } 00792 // konec novinky 00793 00794 00795 for (j=0;j<n;j++){ 00796 ddr[j]=dlambda*v[j]; 00797 ra[j]+=ddr[j]; 00798 fa[j]=fc[j]+(lambda+dlambda)*fp[j]; 00799 } 00800 ddlambda=dlambda; 00801 00802 // computation of internal forces 00803 internal_forces (lcid,fi); 00804 subv(fa, fi, f, n); 00805 norf=sizev(f,n); 00806 norfa = sizev(fa, n); 00807 //if (norfa<1.0e-8){} 00808 //else norf /= norfa; 00809 norf /= norfa; 00810 00811 00812 //if (Mespr==1) fprintf (stdout,"\n ddlambda %e dlambda %e norf %e ierr %e",ddlambda,dlambda,norf,ierr); 00813 if (Mespr==1) fprintf (stdout,"\n increment %ld norv %e norf %e",i,norv,norf); 00814 00815 if (norf<ierr){ 00816 // ****************************************** 00817 // no inner iteration loop is required *** 00818 // ****************************************** 00819 modif++; 00820 00821 if (modif>1){ 00822 // arc length modification 00823 dl*=2.0; 00824 if (dl>dlmax) 00825 dl=dlmax; 00826 modif=0; 00827 if (Mespr==1){ 00828 fprintf (stdout,"\n arc-length must be modified (dl increase) dl=%e",dl); 00829 // fprintf (Out,"\n arc-length must be modified (dl increase) dl=%e",dl); 00830 } 00831 } 00832 lambda+=dlambda; 00833 Mm->updateipval (); 00834 compute_req_val (lcid); 00835 print_step(lcid, i+1, lambda, fa); 00836 print_flush(); 00837 } 00838 else{ 00839 // **************************** 00840 // inner iteration loop **** 00841 // **************************** 00842 stop=0; 00843 for (j=0;j<ini;j++){ 00844 if (Mp->nlman->stmat==tangent_stiff){ 00845 stiffness_matrix (lcid); 00846 Fsd->initialization (Ndofm,Mp->ense,f); 00847 } 00848 00849 fprintf (Out,"\n arc-length vnitrni smycka %ld %ld",i,j); 00850 00851 00852 copyv(f, buf, n); 00853 00854 Fsd->solve_lin_alg_system (u,f); 00855 00856 // back substitution 00857 //Smat->solve_system (u,f); 00858 00859 // backup of the vector ddr 00860 copyv(ddr, buddr, n); 00861 00862 00863 copyv(ddr, f, n); 00864 addv(ddr, u, n); 00865 // coefficient of quadratic equation 00866 quadeqcoeff (ddr,v,n,ddlambda,psi,norfp,dl,a0,a1,a2); 00867 00868 // solution of quadratic equation 00869 numr = solv_polynom_2(a2, a1, a0, l1, l2); 00870 switch (numr) 00871 { 00872 case -1 : 00873 fprintf (stderr,"\n\n infinite number of solution of constrained condition in function arclength"); 00874 break; 00875 case 0 : 00876 fprintf (stderr,"\n\n nonsolvable constrained condition in function arclength"); 00877 break; 00878 case 1 : 00879 dlambda = l1; 00880 break; 00881 default : 00882 break; 00883 } 00884 ss1=0.0; ss2=0.0; 00885 ss3=0.0; ss4=0.0; ss5=0.0; 00886 for (k=0;k<n;k++){ 00887 ss1+=(ddr[k]+l1*v[k])*f[k]; 00888 ss2+=(ddr[k]+l2*v[k])*f[k]; 00889 ss3+=(ddr[k]+l1*v[k])*(ddr[k]+l1*v[k]); 00890 ss4+=(ddr[k]+l2*v[k])*(ddr[k]+l2*v[k]); 00891 ss5+=f[k]*f[k]; 00892 } 00893 00894 if (fabs(l1)>fabs(l2)) dlambda=l2; 00895 else dlambda=l1; 00896 00897 copyv (Fsd->tw,bulm,nlm); 00898 00899 Fsd->add_mult (dlambda); 00900 Fsd->mult_correction (); 00901 00902 // novinka 00903 if (Fsd->nch!=0){ 00904 for (l=0;l<ini;l++){ 00905 00906 fprintf (stdout,"\n cyklime v casti dve %ld",j); 00907 00908 00909 // restoring of Lagrange multipliers 00910 copyv (bulm,Fsd->tw,nlm); 00911 copyv (buddr,ddr,n); 00912 copyv (buf,f,n); 00913 00914 Fsd->solve_lin_alg_system (u,f); 00915 00916 copyv(ddr, f, n); 00917 addv(ddr, u, n); 00918 // coefficient of quadratic equation 00919 quadeqcoeff (ddr,v,n,ddlambda,psi,norfp,dl,a0,a1,a2); 00920 00921 // solution of quadratic equation 00922 numr = solv_polynom_2(a2, a1, a0, l1, l2); 00923 switch (numr) 00924 { 00925 case -1 : 00926 fprintf (stderr,"\n\n infinite number of solution of constrained condition in function arclength"); 00927 break; 00928 case 0 : 00929 fprintf (stderr,"\n\n nonsolvable constrained condition in function arclength"); 00930 break; 00931 case 1 : 00932 dlambda = l1; 00933 break; 00934 default : 00935 break; 00936 } 00937 ss1=0.0; ss2=0.0; 00938 ss3=0.0; ss4=0.0; ss5=0.0; 00939 for (k=0;k<n;k++){ 00940 ss1+=(ddr[k]+l1*v[k])*f[k]; 00941 ss2+=(ddr[k]+l2*v[k])*f[k]; 00942 ss3+=(ddr[k]+l1*v[k])*(ddr[k]+l1*v[k]); 00943 ss4+=(ddr[k]+l2*v[k])*(ddr[k]+l2*v[k]); 00944 ss5+=f[k]*f[k]; 00945 } 00946 00947 if (fabs(l1)>fabs(l2)) dlambda=l2; 00948 else dlambda=l1; 00949 00950 00951 Fsd->add_mult (dlambda); 00952 Fsd->mult_correction (); 00953 00954 if (Fsd->nch==0) 00955 break; 00956 } 00957 } 00958 // konec novinky 00959 00960 00961 00962 00963 for (k=0;k<n;k++){ 00964 ddr[k]+=dlambda*v[k]; 00965 ra[k]+=u[k]+dlambda*v[k]; 00966 fa[k]+=dlambda*fp[k]; 00967 } 00968 ddlambda+=dlambda; 00969 00970 //fprintf (stdout," ddlambda %e",ddlambda); 00971 00972 //fprintf (Out,"\n ddlambda %e dlambda %e",ddlambda,dlambda); 00973 // computation of internal forces 00974 internal_forces (lcid,fi); 00975 subv(fa, fi, f, n); 00976 norf=sizev(f,n); 00977 norfa = sizev(fa, n); 00978 norf /= norfa; 00979 //if (norfa<1.0e-8){} 00980 //else norf /= norfa; 00981 00982 if (Mespr==1){ 00983 fprintf (stdout,"\n norf=%e ierr=%e",norf,ierr); 00984 //fprintf (Out,"\n increment %ld inner loop %ld norf=%e",i,j,norf); 00985 } 00986 00987 if (norf<ierr){ 00988 lambda+=ddlambda; 00989 Mm->updateipval (); 00990 stop=1; 00991 compute_req_val (lcid); 00992 print_step(lcid, i+1, lambda, fa); 00993 print_flush(); 00994 break; 00995 } 00996 } 00997 modif=0; 00998 if (stop==0){ 00999 // modification of the arc length 01000 dl/=2.0; 01001 if (dl<dlmin){ 01002 dl=dlmin; 01003 break; 01004 } 01005 if (Mespr==1){ 01006 fprintf (stdout,"\n arc length must be modified (dl decrease) dl=%e",dl); 01007 //fprintf (Out,"\n arc length must be modified (dl decrease) dl=%e",dl); 01008 } 01009 // restoring of left hand side vector 01010 copyv(r, ra, n); 01011 // restoring of lambda parameter 01012 lambda=blambda; 01013 // restoring of Lagrange multipliers 01014 copyv (blm,Fsd->tw,nlm); 01015 } 01016 } 01017 01018 fprintf (stdout,"\n increment %ld total lambda %e",i,lambda); 01019 01020 if (stop==0) 01021 continue; 01022 01023 } 01024 01025 // ------------------------------------ 01026 // finish of main iteration loop ---- 01027 // ------------------------------------ 01028 01029 01030 fprintf (stdout,"\n\n multiplikatory %lf",Fsd->tw[0]); 01031 01032 if (Mp->nlman->hdbackupal==1) 01033 arclsave (i,n,lambda,dl,ra,fp); 01034 01035 print_close(); 01036 delete [] r; 01037 delete [] fi; 01038 delete [] fa; 01039 delete [] f; 01040 delete [] v; 01041 delete [] u; 01042 delete [] ddr; 01043 01044 } 01045 */ 01046 01047 01048 /** 01049 JK, 26.6.2006 01050 */ 01051 /* 01052 void solve_nonlinear_floating_subdomains_29_6 () 01053 { 01054 long i,j,k,n,lcid,ni,ini,stop,modif,li,numr,nlm; 01055 double a0,a1,a2,l1,l2,dl,dlmax,dlmin,psi,lambda,blambda,dlambda,ddlambda,norf,norfp,norv,norfa,zero,ierr; 01056 double lambdao;//,ss1,ss2,ss3,ss4,ss5; 01057 double *r,*ra,*ddr,*u,*v,*f,*fa,*fp,*fc,*fi,*d,*rhs; 01058 01059 01060 // load case id must be equal to zero 01061 // only one load case can be solved 01062 // one load case may contain several proportional vectors 01063 lcid = 0; 01064 01065 // assembling of stiffness matrix 01066 stiffness_matrix (lcid); 01067 01068 // number of unknown displacements 01069 n = Ndofm; 01070 01071 // auxiliary assigment of pointer 01072 rhs=Lsrs->give_rhs (lcid); 01073 01074 // assembling of right hand side vector (vector of nodal forces) 01075 mefel_right_hand_side (lcid,rhs); 01076 01077 // vector of nodal displacements 01078 d = Lsrs->give_lhs (lcid); 01079 // vector of proportional load 01080 fp = Lsrs->give_rhs (lcid*2); 01081 // vector of constant load 01082 fc = Lsrs->give_rhs (lcid*2+1); 01083 01084 01085 01086 // allocation of object floating subdomains 01087 Fsd = new flsubdom; 01088 01089 // computation of rigid body motions 01090 // list of dependent equations 01091 // determination of number of Lagrange multipliers 01092 Fsd->initialization (Ndofm,Mp->ense,fp); 01093 01094 // number of Lagrange multipliers 01095 nlm = Fsd->nlm; 01096 01097 01098 01099 // ******************************* 01100 // data about arc-length solver 01101 // ******************************* 01102 01103 // maximum number of increments 01104 ni = Mp->nlman->nial; 01105 // maximum number of iterations in one increment 01106 ini = Mp->nlman->niilal; 01107 // computer zero 01108 zero = Mp->zero; 01109 // required error in inner loop 01110 ierr = Mp->nlman->erral; 01111 // length of arc 01112 dl = Mp->nlman->dlal; 01113 // maximum length of arc 01114 dlmax = Mp->nlman->dlmaxal; 01115 // minimum length of arc 01116 dlmin = Mp->nlman->dlminal; 01117 // displacement-loading driving switch 01118 psi = Mp->nlman->psial; 01119 01120 01121 // ********************************* 01122 // allocation of auxiliary arrays 01123 // ********************************* 01124 r = new double [nlm]; 01125 ra = new double [nlm]; 01126 ddr = new double [nlm]; 01127 u = new double [nlm]; 01128 v = new double [nlm]; 01129 f = new double [n]; 01130 fa = new double [n]; 01131 fi = new double [n]; 01132 01133 01134 lambda=0.0; 01135 lambdao=0.0; 01136 dlambda=0.0; 01137 li=0; 01138 modif=0; 01139 01140 01141 01142 01143 01144 // norm of proportionality vector 01145 norfp = ss(fp,fp,n); 01146 01147 01148 01149 // assembling reached load vector 01150 for (j=0;j<n;j++){ 01151 fa[j]=fc[j]+(lambda+dlambda)*fp[j]; 01152 } 01153 if (li) 01154 print_init(-1, "at"); 01155 else 01156 { 01157 print_init(-1, "wt"); 01158 print_step(lcid, li, lambda, fa); 01159 } 01160 01161 01162 01163 // *************************** 01164 // main iteration loop **** 01165 // *************************** 01166 for (i=li;i<ni;i++){ 01167 01168 fprintf (Out,"\n\n arc-length prirustek %ld",i); 01169 01170 01171 stop=1; // termit 01172 // backup of left hand side vector 01173 copyv (ra, r, nlm); 01174 // backup of reached lambda parameter 01175 blambda=lambda; 01176 01177 // backup of the fp, in ldl_sky the right hand side will be destroyed 01178 copyv(fp, f, n); 01179 01180 // solution of K(r).v=F 01181 Fsd->pmcg (f,Out); 01182 copyv (Fsd->w,v,nlm); 01183 01184 // generalized norm of displacement increments 01185 norv = displincr (v,nlm); 01186 01187 // compute new dlambda increment 01188 dlambda = dl/sqrt(norv+psi*psi*norfp); 01189 01190 for (j=0;j<nlm;j++){ 01191 ddr[j]=dlambda*v[j]; 01192 ra[j]+=ddr[j]; 01193 } 01194 for (j=0;j<n;j++){ 01195 fa[j]=fc[j]+(lambda+dlambda)*fp[j]; 01196 } 01197 ddlambda=dlambda; 01198 01199 Fsd->add_mult (dlambda); 01200 Fsd->mult_correction (); 01201 01202 Fsd->nonlinlagrmultdispl (d,fa); 01203 01204 01205 // computation of internal forces 01206 internal_forces (lcid,fi); 01207 subv(fa, fi, f, n); 01208 norf = sizev(f,n); 01209 norfa = sizev(fa, n); 01210 norf /= norfa; 01211 01212 01213 if (Mespr==1) fprintf (stdout,"\n increment %ld norv %e norf %e",i,norv,norf); 01214 01215 if (norf<ierr){ 01216 // ****************************************** 01217 // no inner iteration loop is required *** 01218 // ****************************************** 01219 modif++; 01220 01221 if (modif>1){ 01222 // arc length modification 01223 dl*=2.0; 01224 if (dl>dlmax) 01225 dl=dlmax; 01226 modif=0; 01227 if (Mespr==1){ 01228 fprintf (stdout,"\n arc-length must be modified (dl increase) dl=%e",dl); 01229 // fprintf (Out,"\n arc-length must be modified (dl increase) dl=%e",dl); 01230 } 01231 } 01232 lambda+=dlambda; 01233 Mm->updateipval (); 01234 compute_req_val (lcid); 01235 print_step(lcid, i+1, lambda, fa); 01236 print_flush(); 01237 } 01238 else{ 01239 // **************************** 01240 // inner iteration loop **** 01241 // **************************** 01242 stop=0; 01243 for (j=0;j<ini;j++){ 01244 01245 fprintf (Out,"\n arc-length vnitrni smycka %ld %ld",i,j); 01246 01247 // solution of K(r).v=F 01248 Fsd->pmcg (f,Out); 01249 copyv (Fsd->w,u,nlm); 01250 01251 01252 //copyv(ddr, f, nlm); 01253 addv(ddr, u, nlm); 01254 // coefficient of quadratic equation 01255 //quadeqcoeff (ddr,v,n,ddlambda,psi,norfp,dl,a0,a1,a2); 01256 quadeqcoeff (ddr,v,nlm,ddlambda,psi,norfp,dl,a0,a1,a2); 01257 01258 fprintf (Out,"\n\n\n Kontrola kvadraticke rovnice"); 01259 fprintf (Out,"\n norfp %15.10le",norfp); 01260 fprintf (Out,"\n norv %15.10le",norv); 01261 fprintf (Out,"\n (ddr,v) %15.10le",ss(ddr,v,n)); 01262 fprintf (Out,"\n (ddr,ddr) %15.10le",ss(ddr,ddr,n)); 01263 fprintf (Out,"\n dl %15.10le",dl); 01264 fprintf (Out,"\n ddlambda %15.10le",ddlambda); 01265 fprintf (Out,"\n psi %15.10le",psi); 01266 fprintf (Out,"\n a2 %15.10le",a2); 01267 fprintf (Out,"\n a1 %15.10le",a1); 01268 fprintf (Out,"\n a0 %15.10le",a0); 01269 fprintf (Out,"\n discrim %15.10le",a1*a1-4.0*a2*a0); 01270 fprintf (Out,"\n\n"); 01271 01272 01273 // solution of quadratic equation 01274 numr = solv_polynom_2(a2, a1, a0, l1, l2); 01275 switch (numr) 01276 { 01277 case -1 : 01278 fprintf (stderr,"\n\n infinite number of solution of constrained condition in function arclength"); 01279 break; 01280 case 0 : 01281 fprintf (stderr,"\n\n nonsolvable constrained condition in function arclength"); 01282 break; 01283 case 1 : 01284 dlambda = l1; 01285 break; 01286 default : 01287 break; 01288 } 01289 01290 if (fabs(l1)>fabs(l2)) dlambda=l2; 01291 else dlambda=l1; 01292 01293 for (k=0;k<nlm;k++){ 01294 ddr[k]+=dlambda*v[k]; 01295 ra[k]+=u[k]+dlambda*v[k]; 01296 } 01297 for (k=0;k<n;k++){ 01298 fa[k]+=dlambda*fp[k]; 01299 } 01300 ddlambda+=dlambda; 01301 01302 Fsd->add_mult (dlambda); 01303 Fsd->mult_correction (); 01304 01305 Fsd->nonlinlagrmultdispl (d,fa); 01306 01307 // computation of internal forces 01308 internal_forces (lcid,fi); 01309 subv(fa, fi, f, n); 01310 norf=sizev(f,n); 01311 norfa = sizev(fa, n); 01312 norf /= norfa; 01313 01314 if (Mespr==1){ 01315 fprintf (stdout,"\n norf=%e ierr=%e",norf,ierr); 01316 //fprintf (Out,"\n increment %ld inner loop %ld norf=%e",i,j,norf); 01317 } 01318 01319 if (norf<ierr){ 01320 lambda+=ddlambda; 01321 Mm->updateipval (); 01322 stop=1; 01323 compute_req_val (lcid); 01324 print_step(lcid, i+1, lambda, fa); 01325 print_flush(); 01326 break; 01327 } 01328 } 01329 modif=0; 01330 if (stop==0){ 01331 // modification of the arc length 01332 dl/=2.0; 01333 if (dl<dlmin){ 01334 dl=dlmin; 01335 break; 01336 } 01337 if (Mespr==1){ 01338 fprintf (stdout,"\n arc length must be modified (dl decrease) dl=%e",dl); 01339 //fprintf (Out,"\n arc length must be modified (dl decrease) dl=%e",dl); 01340 } 01341 // restoring of left hand side vector 01342 copyv(r, ra, nlm); 01343 // restoring of lambda parameter 01344 lambda=blambda; 01345 } 01346 } 01347 01348 fprintf (stdout,"\n increment %ld total lambda %e",i,lambda); 01349 01350 if (stop==0) 01351 continue; 01352 01353 } 01354 01355 // ------------------------------------ 01356 // finish of main iteration loop ---- 01357 // ------------------------------------ 01358 01359 01360 print_close(); 01361 delete [] r; 01362 delete [] fi; 01363 delete [] fa; 01364 delete [] f; 01365 delete [] v; 01366 delete [] u; 01367 delete [] ddr; 01368 01369 } 01370 */ 01371 01372 /* 01373 void solve_nonlinear_floating_subdomains () 01374 { 01375 long i,j,k,n,lcid,nlm,ni,ini,stop,modif,li,newmesh, numr; 01376 double a0,a1,a2,l1,l2,dl,dlmax,dlmin,psi,lambda,blambda,dlambda,ddlambda,norf,norfp,norv,norfa,zero,ierr; 01377 double lambdao,ss1,ss2,ss3,ss4,ss5; 01378 double *r,*ra,*ddr,*u,*v,*f,*fa,*fp,*fc,*fi,*rhs; 01379 01380 01381 // load case id must be equal to zero 01382 // only one load case can be solved 01383 // one load case may contain several proportional vectors 01384 lcid = 0; 01385 01386 // assembling of stiffness matrix 01387 stiffness_matrix (lcid); 01388 01389 // number of unknown displacements 01390 n = Ndofm; 01391 01392 // auxiliary assigment of pointer 01393 rhs=Lsrs->give_rhs (lcid); 01394 01395 // assembling of right hand side vector (vector of nodal forces) 01396 mefel_right_hand_side (lcid,rhs); 01397 01398 // vector of nodal displacements 01399 ra = Lsrs->give_lhs (lcid); 01400 // vector of proportional load 01401 fp = Lsrs->give_rhs (lcid*2); 01402 // vector of constant load 01403 fc = Lsrs->give_rhs (lcid*2+1); 01404 01405 01406 01407 // allocation of object floating subdomains 01408 Fsd = new flsubdom; 01409 01410 // computation of rigid body motions 01411 // list of dependent equations 01412 // determination of number of Lagrange multipliers 01413 Fsd->initialization (Ndofm,Mp->ense,fp); 01414 01415 // number of Lagrange multipliers 01416 nlm = Fsd->nlm; 01417 01418 01419 01420 // ******************************* 01421 // data about arc-length solver 01422 // ******************************* 01423 01424 // maximum number of increments 01425 ni = Mp->nlman->nial; 01426 // maximum number of iterations in one increment 01427 ini = Mp->nlman->niilal; 01428 // computer zero 01429 zero = Mp->zero; 01430 // required error in inner loop 01431 ierr = Mp->nlman->erral; 01432 // length of arc 01433 dl = Mp->nlman->dlal; 01434 // maximum length of arc 01435 dlmax = Mp->nlman->dlmaxal; 01436 // minimum length of arc 01437 dlmin = Mp->nlman->dlminal; 01438 // displacement-loading driving switch 01439 psi = Mp->nlman->psial; 01440 01441 01442 01443 // allocation of auxiliary arrays 01444 r = new double [n]; 01445 ddr = new double [n]; 01446 u = new double [n]; 01447 v = new double [n]; 01448 f = new double [n]; 01449 fa = new double [n]; 01450 fi = new double [n]; 01451 01452 01453 dlambda=0.0; 01454 lambda=0.0; 01455 lambdao=0.0; 01456 li=0; 01457 01458 01459 // norm of proportionality vector 01460 norfp = ss(fp,fp,n); 01461 modif=0; 01462 01463 01464 // assembling reached load vector 01465 for (j=0;j<n;j++){ 01466 fa[j]=fc[j]+(lambda+dlambda)*fp[j]; 01467 } 01468 if (li) 01469 print_init(-1, "at"); 01470 else 01471 { 01472 print_init(-1, "wt"); 01473 print_step(lcid, li, lambda, fa); 01474 } 01475 01476 01477 01478 // *************************** 01479 // main iteration loop **** 01480 // *************************** 01481 for (i=li;i<ni;i++){ 01482 01483 fprintf (Out,"\n\n arc-length prirustek %ld",i); 01484 01485 01486 stop=1; // termit 01487 // backup of left hand side vector 01488 copyv (ra, r, n); 01489 // backup of reached lambda parameter 01490 blambda=lambda; 01491 01492 01493 // backup of the fp, in ldl_sky the right hand side will be destroyed 01494 copyv(fp, f, n); 01495 01496 01497 stiffness_matrix (lcid); 01498 Fsd->factorize (); 01499 01500 01501 // solution of K(r).v=F 01502 Fsd->pmcg (f,Out); 01503 01504 Fsd->lagrmultdispl (v,f); 01505 01506 01507 // generalized norm of displacement increments 01508 norv = displincr (v,n); 01509 01510 // compute new dlambda increment 01511 dlambda = dl/sqrt(norv+psi*psi*norfp); 01512 01513 //if (norv+psi*psi*norfp < 1.0e-8) dlambda=0.0; 01514 //else dlambda = dl/sqrt(norv+psi*psi*norfp); 01515 01516 for (j=0;j<n;j++){ 01517 ddr[j]=dlambda*v[j]; 01518 ra[j]+=ddr[j]; 01519 fa[j]=fc[j]+(lambda+dlambda)*fp[j]; 01520 } 01521 ddlambda=dlambda; 01522 01523 Fsd->add_mult (dlambda); 01524 Fsd->mult_correction (); 01525 01526 // computation of internal forces 01527 internal_forces (lcid,fi); 01528 subv(fa, fi, f, n); 01529 norf=sizev(f,n); 01530 norfa = sizev(fa, n); 01531 //if (norfa<1.0e-8){} 01532 //else norf /= norfa; 01533 norf /= norfa; 01534 01535 01536 //if (Mespr==1) fprintf (stdout,"\n ddlambda %e dlambda %e norf %e ierr %e",ddlambda,dlambda,norf,ierr); 01537 if (Mespr==1) fprintf (stdout,"\n increment %ld norv %e norf %e",i,norv,norf); 01538 01539 if (norf<ierr){ 01540 // ****************************************** 01541 // no inner iteration loop is required *** 01542 // ****************************************** 01543 modif++; 01544 01545 if (modif>1){ 01546 // arc length modification 01547 dl*=2.0; 01548 if (dl>dlmax) 01549 dl=dlmax; 01550 modif=0; 01551 if (Mespr==1){ 01552 fprintf (stdout,"\n arc-length must be modified (dl increase) dl=%e",dl); 01553 // fprintf (Out,"\n arc-length must be modified (dl increase) dl=%e",dl); 01554 } 01555 } 01556 lambda+=dlambda; 01557 Mm->updateipval (); 01558 compute_req_val (lcid); 01559 print_step(lcid, i+1, lambda, fa); 01560 print_flush(); 01561 } 01562 else{ 01563 // **************************** 01564 // inner iteration loop **** 01565 // **************************** 01566 stop=0; 01567 for (j=0;j<ini;j++){ 01568 01569 01570 fprintf (Out,"\n arc-length vnitrni smycka %ld %ld",i,j); 01571 fprintf (stdout,"\n arc-length vnitrni smycka %ld %ld",i,j); 01572 01573 01574 stiffness_matrix (lcid); 01575 Fsd->factorize (); 01576 01577 01578 // back substitution 01579 Fsd->pmcg (f,Out); 01580 Fsd->lagrmultdispl (u,f); 01581 01582 copyv(ddr, f, n); 01583 addv(ddr, u, n); 01584 // coefficient of quadratic equation 01585 quadeqcoeff (ddr,v,n,ddlambda,psi,norfp,dl,a0,a1,a2); 01586 01587 // solution of quadratic equation 01588 numr = solv_polynom_2(a2, a1, a0, l1, l2); 01589 switch (numr) 01590 { 01591 case -1 : 01592 fprintf (stderr,"\n\n infinite number of solution of constrained condition in function arclength"); 01593 break; 01594 case 0 : 01595 fprintf (stderr,"\n\n nonsolvable constrained condition in function arclength"); 01596 break; 01597 case 1 : 01598 dlambda = l1; 01599 break; 01600 default : 01601 break; 01602 } 01603 01604 if (fabs(l1)>fabs(l2)) dlambda=l2; 01605 else dlambda=l1; 01606 01607 //fprintf (stdout,"\n increment %ld inner loop %ld x1=%e x2=%e dlambda %e",i,j,l1,l2,dlambda); 01608 //fprintf (stdout,"\n increment %ld inner loop %ld x1=%e x2=%e ddlambda %e",i,j,l1,l2,ddlambda); 01609 01610 for (k=0;k<n;k++){ 01611 ddr[k]+=dlambda*v[k]; 01612 ra[k]+=u[k]+dlambda*v[k]; 01613 fa[k]+=dlambda*fp[k]; 01614 } 01615 ddlambda+=dlambda; 01616 01617 Fsd->add_mult (dlambda); 01618 Fsd->mult_correction (); 01619 01620 01621 // computation of internal forces 01622 internal_forces (lcid,fi); 01623 subv(fa, fi, f, n); 01624 norf=sizev(f,n); 01625 norfa = sizev(fa, n); 01626 norf /= norfa; 01627 01628 if (Mespr==1){ 01629 fprintf (stdout,"\n norf=%e ierr=%e",norf,ierr); 01630 //fprintf (Out,"\n increment %ld inner loop %ld norf=%e",i,j,norf); 01631 } 01632 01633 if (norf<ierr){ 01634 lambda+=ddlambda; 01635 Mm->updateipval (); 01636 stop=1; 01637 compute_req_val (lcid); 01638 print_step(lcid, i+1, lambda, fa); 01639 print_flush(); 01640 break; 01641 } 01642 } 01643 modif=0; 01644 if (stop==0){ 01645 // modification of the arc length 01646 dl/=2.0; 01647 if (dl<dlmin){ 01648 dl=dlmin; 01649 break; 01650 } 01651 if (Mespr==1){ 01652 fprintf (stdout,"\n arc length must be modified (dl decrease) dl=%e",dl); 01653 //fprintf (Out,"\n arc length must be modified (dl decrease) dl=%e",dl); 01654 } 01655 // restoring of left hand side vector 01656 copyv(r, ra, n); 01657 // restoring of lambda parameter 01658 lambda=blambda; 01659 } 01660 } 01661 01662 fprintf (stdout,"\n increment %ld total lambda %e",i,lambda); 01663 01664 if (stop==0) 01665 continue; 01666 01667 } 01668 01669 // ------------------------------------ 01670 // finish of main iteration loop ---- 01671 // ------------------------------------ 01672 01673 01674 if (Mp->nlman->hdbackupal==1) 01675 arclsave (i,n,lambda,dl,ra,fp); 01676 01677 print_close(); 01678 delete [] r; 01679 delete [] fi; 01680 delete [] fa; 01681 delete [] f; 01682 delete [] v; 01683 delete [] u; 01684 delete [] ddr; 01685 01686 } 01687 */ 01688 01689 01690 /** 01691 function solves system of nonlinear algebraic equation 01692 by arc-length method (continuation method) 01693 only one right hand side vector is supported with respect 01694 to nonlinearity and absence of superposition 01695 01696 d stands for delta 01697 dd stands for capital delta 01698 01699 fc - nonproportional %vector 01700 fp - proportional %vector 01701 n - order of the system 01702 01703 @param lcid - load case id 01704 01705 25.7.2001 01706 */ 01707 /* 01708 void solve_nonlinear_floating_subdomains () 01709 { 01710 long i,j,k,n,nlm,ni,ini,li,numr,stop,modif,lcid; 01711 double dl,dlmax,dlmin,psi,ierr,zero,lambda,dlambda,ddlambda,blambda,norf,norfa; 01712 double a0,a1,a2,l1,l2,norv; 01713 double *r,*ddr,*fc,*fp,*f,*dv,*du,*dnu,*rhs; 01714 double *br,*bmu,*av,*av1,*av2,*fa,*fi,*aa1,*aa2,*a,*b,*c; 01715 01716 01717 //long i,j,k,n,ni,ini,stop,modif,li,newmesh, numr; 01718 //double a0,a1,a2,l1,l2,dl,dlmax,dlmin,psi,lambda,blambda,dlambda,ddlambda,norf,norfp,norv,norfa,zero,ierr; 01719 //double lambdao,ss1,ss2,ss3,ss4,ss5; 01720 //double *r,*ra,*ddr,*u,*v,*f,*fa,*fp,*fc,*fi; 01721 01722 01723 01724 // number of rows of the matrix 01725 n = Ndofm; 01726 // maximum number of increments 01727 ni = Mp->nlman->nial; 01728 // maximum number of iterations in one increment 01729 ini = Mp->nlman->niilal; 01730 // computer zero 01731 zero = Mp->zero; 01732 // required error in inner loop 01733 ierr = Mp->nlman->erral; 01734 // length of arc 01735 dl = Mp->nlman->dlal; 01736 // maximum length of arc 01737 dlmax = Mp->nlman->dlmaxal; 01738 // minimum length of arc 01739 dlmin = Mp->nlman->dlminal; 01740 // displacement-loading driving switch 01741 psi = Mp->nlman->psial; 01742 01743 // load case id must be equal to zero 01744 // only one load case can be solved 01745 // one load case may contain several proportional vectors 01746 lcid = 0; 01747 01748 // assembling of stiffness matrix 01749 stiffness_matrix (lcid); 01750 01751 01752 // auxiliary assigment of pointer 01753 rhs=Lsrs->give_rhs (lcid); 01754 01755 // assembling of right hand side vector (vector of nodal forces) 01756 mefel_right_hand_side (lcid,rhs); 01757 01758 // vector of nodal displacements 01759 r = Lsrs->give_lhs (lcid); 01760 // vector of proportional load 01761 fp = Lsrs->give_rhs (lcid*2); 01762 // vector of constant load 01763 fc = Lsrs->give_rhs (lcid*2+1); 01764 01765 01766 01767 // allocation of object floating subdomains 01768 Fsd = new flsubdom; 01769 01770 // computation of rigid body motions 01771 // list of dependent equations 01772 // determination of number of Lagrange multipliers 01773 Fsd->initialization (Ndofm,Mp->ense,fp); 01774 01775 // number of Lagrange multipliers 01776 nlm = Fsd->nlm; 01777 01778 01779 // allocation of auxiliary arrays 01780 // array for backup of nodal displacements 01781 br = new double [n]; 01782 // array for backup of Lagrange multipliers 01783 bmu = new double [nlm]; 01784 01785 du = new double [n]; 01786 dv = new double [n]; 01787 f = new double [n]; 01788 fa = new double [n]; 01789 fi = new double [n]; 01790 av = new double [n]; 01791 av1 = new double [n]; 01792 av2 = new double [n]; 01793 aa1 = new double [n]; 01794 aa2 = new double [n]; 01795 01796 dnu = new double [nlm]; 01797 01798 a = new double [n]; 01799 b = new double [n]; 01800 c = new double [n]; 01801 01802 01803 // proportional factor 01804 lambda = 0.0; 01805 dlambda = 0.0; 01806 ddlambda = 0.0; 01807 01808 // attained level in arclength method 01809 // it is 0 by default 01810 // otherwise it is equal to performed number of steps in previous computations 01811 li=0; 01812 01813 modif=0; 01814 01815 01816 01817 // initiation of load vector 01818 for (j=0;j<n;j++){ 01819 fa[j]=fc[j]+lambda*fp[j]; 01820 } 01821 01822 01823 if (li) 01824 print_init(-1, "at"); 01825 else 01826 { 01827 print_init(-1, "wt"); 01828 print_step(lcid, li, lambda, fa); 01829 } 01830 01831 01832 // *************************** 01833 // main iteration loop **** 01834 // *************************** 01835 for (i=li;i<ni;i++){ 01836 01837 fprintf (Out,"\n\n arc-length increment %ld",i); 01838 fprintf (stdout,"\n\n arc-length increment %ld",i); 01839 01840 01841 //stop=1; // termit 01842 01843 01844 01845 // backup of nodal displacements 01846 copyv (r, br, n); 01847 // backup of Lagrange multipliers 01848 copyv (Fsd->muo,bmu,nlm); 01849 // backup of reached lambda parameter 01850 blambda=lambda; 01851 01852 01853 01854 01855 // assembling of tangent stiffness matrix 01856 if ((Mp->nlman->stmat==tangent_stiff) || (i == li)){ 01857 stiffness_matrix (lcid); 01858 Fsd->factorize (); 01859 } 01860 01861 // solution of K(r).dv = fp-B.dnu 01862 Fsd->pmcg (fp,Out); 01863 01864 for (k=0;k<nlm;k++){ 01865 fprintf (Out,"\n krok %ld lambda %ld %lf",i,k,Fsd->w[k]); 01866 } 01867 01868 01869 Fsd->lagrmultdispl (dv,fp); 01870 copyv (Fsd->w,dnu,nlm); 01871 01872 01873 // generalized norm of displacement increments 01874 norv = displincr (dv,n); 01875 01876 Gtm->flsub.coarse_local (dnu,av); 01877 01878 norf=0.0; 01879 for (k=0;k<n;k++){ 01880 norf+=(fp[k]-av[k])*(fp[k]-av[k]); 01881 } 01882 01883 // compute new dlambda increment 01884 dlambda = dl/sqrt(norv+psi*psi*norf); 01885 01886 01887 for (k=0;k<n;k++){ 01888 ddr[k]=dlambda*dv[k]; 01889 r[k]+=ddr[k]; 01890 fa[k]=fc[k]+(lambda+dlambda)*fp[k]; 01891 } 01892 ddlambda=dlambda; 01893 for (k=0;k<nlm;k++){ 01894 Fsd->ddmu[k]=dlambda*dnu[k]; 01895 Fsd->mu[k]+=Fsd->ddmu[k]; 01896 } 01897 Fsd->mult_correction (); 01898 01899 01900 // computation of internal forces 01901 internal_forces (lcid,fi); 01902 subv(fa, fi, f, n); 01903 norf=sizev(f,n); 01904 norfa = sizev(fa, n); 01905 norf /= norfa; 01906 01907 01908 //if (Mespr==1) fprintf (stdout,"\n ddlambda %e dlambda %e norf %e ierr %e",ddlambda,dlambda,norf,ierr); 01909 if (Mespr==1) fprintf (stdout,"\n increment %ld norv %e norf %e",i,norv,norf); 01910 01911 if (norf<ierr){ 01912 // ****************************************** 01913 // no inner iteration loop is required *** 01914 // ****************************************** 01915 modif++; 01916 01917 if (modif>1){ 01918 // arc length modification 01919 dl*=2.0; 01920 if (dl>dlmax) 01921 dl=dlmax; 01922 modif=0; 01923 if (Mespr==1){ 01924 fprintf (stdout,"\n arc-length must be modified (dl increase) dl=%e",dl); 01925 // fprintf (Out,"\n arc-length must be modified (dl increase) dl=%e",dl); 01926 } 01927 } 01928 lambda+=dlambda; 01929 01930 for (k=0;k<nlm;k++){ 01931 Fsd->muo[k]+=Fsd->ddmu[k]; 01932 } 01933 01934 Mm->updateipval (); 01935 compute_req_val (lcid); 01936 print_step(lcid, i+1, lambda, fa); 01937 print_flush(); 01938 } 01939 else{ 01940 // **************************** 01941 // inner iteration loop **** 01942 // **************************** 01943 stop=0; 01944 for (j=0;j<ini;j++){ 01945 if (Mp->nlman->stmat==tangent_stiff) 01946 stiffness_matrix (lcid); 01947 01948 Fsd->factorize (); 01949 01950 // back substitution 01951 Fsd->pmcg (f,Out); 01952 Fsd->lagrmultdispl (du,f); 01953 01954 for (k=0;k<nlm;k++){ 01955 fprintf (Out,"\n krok %ld iterace %ld lambda %ld %lf",i,j,k,Fsd->w[k]); 01956 } 01957 01958 fprintf (Out,"\n arc-length vnitrni smycka %ld %ld",i,j); 01959 01960 01961 Gtm->flsub.coarse_local (Fsd->ddmu,av1); 01962 Gtm->flsub.coarse_local (Fsd->w,av2); 01963 01964 for (k=0;k<n;k++){ 01965 a[k]=ddr[k]+du[k]; 01966 } 01967 01968 for (k=0;k<n;k++){ 01969 b[k]=ddlambda*fp[k]-av1[k]-av2[k]; 01970 } 01971 01972 01973 Gtm->flsub.coarse_local (dnu,av); 01974 for (k=0;k<n;k++){ 01975 c[k]=fp[k]-av[k]; 01976 } 01977 01978 a2=0.0; a1=0.0; a0=0.0; 01979 for (k=0;k<n;k++){ 01980 a2+=c[k]*c[k]*psi*psi+dv[k]*dv[k]; 01981 a1+=2.0*b[k]*c[k]*psi*psi+2.0*dv[k]*a[k]; 01982 a0+=a[k]*a[k]+b[k]*b[k]*psi*psi; 01983 } 01984 a0-=dl*dl; 01985 01986 // solution of quadratic equation 01987 numr = solv_polynom_2(a2, a1, a0, l1, l2); 01988 switch (numr) 01989 { 01990 case -1 : 01991 fprintf (stderr,"\n\n infinite number of solution of constrained condition in function arclength"); 01992 break; 01993 case 0 : 01994 fprintf (stderr,"\n\n nonsolvable constrained condition in function arclength"); 01995 break; 01996 case 1 : 01997 dlambda = l1; 01998 break; 01999 default : 02000 break; 02001 } 02002 02003 if (fabs(l1)>fabs(l2)) dlambda=l2; 02004 else dlambda=l1; 02005 02006 02007 02008 02009 02010 for (k=0;k<n;k++){ 02011 ddr[k]+=du[k]+dlambda*dv[k]; 02012 r[k]+=du[k]+dlambda*dv[k]; 02013 fa[k]+=dlambda*fp[k]; 02014 } 02015 ddlambda+=dlambda; 02016 for (k=0;k<nlm;k++){ 02017 Fsd->ddmu[k]+=Fsd->w[k]+dlambda*dnu[k]; 02018 Fsd->mu[k]+=Fsd->w[k]+dlambda*dnu[k]; 02019 } 02020 Fsd->mult_correction (); 02021 02022 //fprintf (stdout," ddlambda %e",ddlambda); 02023 02024 //fprintf (Out,"\n ddlambda %e dlambda %e",ddlambda,dlambda); 02025 // computation of internal forces 02026 internal_forces (lcid,fi); 02027 subv(fa, fi, f, n); 02028 norf=sizev(f,n); 02029 norfa = sizev(fa, n); 02030 norf /= norfa; 02031 02032 if (Mespr==1){ 02033 fprintf (stdout,"\n norf=%e ierr=%e",norf,ierr); 02034 //fprintf (Out,"\n increment %ld inner loop %ld norf=%e",i,j,norf); 02035 } 02036 02037 if (norf<ierr){ 02038 lambda+=ddlambda; 02039 02040 for (k=0;k<nlm;k++){ 02041 Fsd->muo[k]+=Fsd->ddmu[k]; 02042 } 02043 02044 Mm->updateipval (); 02045 stop=1; 02046 compute_req_val (lcid); 02047 print_step(lcid, i+1, lambda, fa); 02048 print_flush(); 02049 break; 02050 } 02051 } 02052 modif=0; 02053 if (stop==0){ 02054 // modification of the arc length 02055 dl/=2.0; 02056 if (dl<dlmin){ 02057 dl=dlmin; 02058 break; 02059 } 02060 if (Mespr==1){ 02061 fprintf (stdout,"\n arc length must be modified (dl decrease) dl=%e",dl); 02062 //fprintf (Out,"\n arc length must be modified (dl decrease) dl=%e",dl); 02063 } 02064 // restoring of nodal displacements 02065 copyv(br, r, n); 02066 // restoring of Lagrange multipliers 02067 copyv(bmu,Fsd->muo,nlm); 02068 // restoring of lambda parameter 02069 lambda=blambda; 02070 } 02071 } 02072 02073 fprintf (stdout,"\n increment %ld total lambda %e",i,lambda); 02074 02075 if (stop==0) 02076 continue; 02077 02078 } 02079 02080 // ------------------------------------ 02081 // finish of main iteration loop ---- 02082 // ------------------------------------ 02083 02084 02085 02086 02087 //if (Mp->nlman->hdbackupal==1) 02088 //arclsave (i,n,lambda,dl,ra,fp); 02089 02090 print_close(); 02091 02092 02093 delete [] br; 02094 delete [] bmu; 02095 delete [] dv; 02096 delete [] f; 02097 delete [] fa; 02098 delete [] fi; 02099 delete [] av; 02100 delete [] av1; 02101 delete [] av2; 02102 delete [] aa1; 02103 delete [] aa2; 02104 delete [] dnu; 02105 02106 } 02107 */ 02108 02109 02110 02111 02112 /** 02113 function solves system of nonlinear algebraic equations 02114 load is applied by prescribed displacements 02115 02116 JK, 19.1.2007 02117 */ 02118 /* 02119 void solve_nonlinear_floating_subdomains_19_1 () 02120 { 02121 long i,j,ni,ini,lcid; 02122 double lambda,dlambda,mindlambda,err,norfi; 02123 double *f,*lhs,*rhs,*arhs,*fi; 02124 02125 // only one load case is assumed 02126 lcid=0; 02127 // number of iteration steps/increments 02128 ni=Mp->nlman->ninr; 02129 // maximum number of iterations in one increment 02130 ini = Mp->nlman->niilnr; 02131 // total magnitude of the right hand side 02132 lambda=0.0; 02133 // magnitude of the step increment 02134 dlambda=Mp->nlman->incrnr; 02135 // minimum magnitude of the step increment 02136 mindlambda=Mp->nlman->minincrnr; 02137 // required error 02138 err=Mp->nlman->errnr; 02139 02140 // 02141 // fi = new double [?] 02142 02143 Fsd = new flsubdom; 02144 02145 // nodal displacements 02146 lhs=Lsrs->give_lhs (lcid); 02147 // proportional vector of the prescribed forces 02148 rhs=Lsrs->give_rhs (lcid); 02149 // auxiliary array 02150 f = new double [Ndofm]; 02151 // auxiliary vector 02152 arhs = new double [Ndofm]; 02153 02154 02155 // proportional vector of the prescribed forces 02156 mefel_right_hand_side (lcid,rhs); 02157 02158 // stiffness matrix 02159 stiffness_matrix (lcid); 02160 02161 Fsd->initialization (Ndofm,Mp->ense,rhs); 02162 02163 // loop over iteration steps/increments 02164 for (i=0;i<ni;i++){ 02165 02166 // increment of the right hand side 02167 for (j=0;j<Ndofm;j++){ 02168 f[j]=dlambda*rhs[j]; 02169 } 02170 02171 // solution of the system of algebraic equations 02172 // increment of the Lagrange multipliers are obtained 02173 Fsd->pmcg (f,Out); 02174 02175 // update of Lagrange multipliers 02176 Fsd->add_mult (1.0); 02177 02178 // correction of the Lagrange multipliers 02179 Fsd->mult_correction (); 02180 02181 02182 for (j=0;j<Ndofm;j++){ 02183 arhs[j]=lambda*rhs[j]; 02184 } 02185 02186 // computation of displacements from Lagrange multipliers 02187 Fsd->lagrmultdispl (lhs,arhs); 02188 02189 02190 // actual internal forces 02191 internal_forces (lcid,fi); 02192 02193 // norm of the vector of internal forces 02194 norfi = ss (fi,fi,Ndofm); 02195 02196 if (norfi<err){ 02197 // no inner loop is required 02198 } 02199 else{ 02200 // inner loop is required 02201 for (j=0;j<ini;j++){ 02202 02203 Fsd->solve_lin_alg_system (lhs,fi); 02204 02205 02206 internal_forces (lcid,fi); 02207 02208 } 02209 02210 } 02211 02212 02213 print_init(-1, "wt"); 02214 for (i=0;i<Lsrs->nlc;i++){ 02215 //compute_ipstrains (i); 02216 //compute_ipstresses (i); 02217 //compute_req_val (i); 02218 print_step(i, 0, 0.0, NULL); 02219 } 02220 02221 02222 } 02223 02224 print_close(); 02225 02226 02227 } 02228 */