#include "epsolver.h"#include "arclength.h"#include "global.h"#include "globmat.h"#include "mechprint.h"#include "loadcase.h"#include "dloadcase.h"#include "nssolver.h"#include "springel.h"#include "intpoints.h"#include "vector.h"#include "matrix.h"#include "mathem.h"#include "gmatrix.h"#include "tensor.h"#include "vecttens.h"#include <math.h>#include <string.h>Go to the source code of this file.
Functions | |
| double | arclengthrv (long lcid, double rlambda, double rlerr) |
| the arc-length method which stops iteration on the given value of the lambda parameter | |
| void | epressure_solver (long lcid) |
| function controling earth pressure solvers | |
| void | femplast_epressure (long lcid) |
| earth pressure solver with temporary supports for plasticity models | |
| void | general_epressure (long lcid) |
| earth pressure solver with load iteration | |
| void | general_epressure_varsup (long lcid) |
| earth pressure solver with temporary supports | |
| double | newton_raphsonep (long lcid, double *pfp) |
| the Newton-Raphsonmethop modified for starting from a reached state | |
| void | restore_displ (long lcid, double **bckr) |
| extraction of the reached displacement form the backup vector | |
| void | solve_epressure () |
| main computation function for earth pressure problem type | |
| double arclengthrv | ( | long | lcid, | |
| double | rlambda, | |||
| double | rlerr | |||
| ) |
the arc-length method which stops iteration on the given value of the lambda parameter
The function solves system of nonlinear algebraic equations by the arc-length method. This method was modified to stop iteration on the given value rlambda of the lambda parameter. Rlambda value is reached with required error rlerr. Only one right hand side vector is supported with respect to nonlinearity and absence of superposition.
| lcid | - load case id | |
| rlambda | - required value of lambda parameter | |
| rlerr | - required error between lambda and rlambda |
Cretated by Tomas Koudelka
Definition at line 465 of file epsolver.cpp.
References addv(), arclopen(), arclsave(), copyv(), displincr(), nonlinman::dlal, nonlinman::dlmaxal, nonlinman::dlminal, nonlinman::erral, f, file, probdesc::filename, mechtop::give_code_numbers(), mechtop::give_dof(), mechtop::give_elem_type(), lhsrhs::give_lhs(), lhsrhs::give_rhs(), Gtm, nonlinman::hdbackupal, internal_forces(), loadincr(), Lsrs, Mespr, Mm, mode, Mp, Mt, Ndofm, mechtop::ne, nonlinman::nial, nonlinman::niilal, probdesc::nlman, mechtop::nn, normv(), Out, probdesc::path, mechmat::plast, print_err(), print_flush(), print_step(), nonlinman::psial, quadeqcoeff(), springel::res_stiffness_matrix(), Smat, solv_polynom_2(), slesolv::solve_system(), Spring, spring_1, probdesc::ssle, stiffness_matrix(), subv(), mechmat::updateipval(), and probdesc::zero.
| void epressure_solver | ( | long | lcid | ) |
function controling earth pressure solvers
The function calls appropriate solver for the earth pressure problem.
| lcid | - load case id. |
Created by Tomas Koudelka
Definition at line 67 of file epsolver.cpp.
References epplast_sol, femplast_epressure(), general_epressure(), general_epressure_varsup(), gep_sol, gepvarsup_sol, Mp, print_err(), and probdesc::tepsol.
Referenced by solve_epressure().
| void femplast_epressure | ( | long | lcid | ) |
earth pressure solver with temporary supports for plasticity models
The function solves erath pressure problem using FEM plasticity models. Iteration of rigid spring support is used. Spring supports are incrementally removed during each itreation step. Removed supports are expected to be ordered as last elements, each step is removed the last element and number of elements is decremented. So the user should order removed spring elements so the first removed spring element is the last element (element with the greatest number) and so on. Newton-Raphson method is used for given plasticity problem solution. This solver supposes the static load equal to dead weight of excaved soil is applied to the bottom of excaved part of the soil body and this static load is decremented in the each step when the spring support is removed.
| lcid | - load case id. |
Created by Tomas Koudelka
Definition at line 307 of file epsolver.cpp.
References ivector::a, addv(), allocv(), mechmat::computenlstresses(), copyv(), destrv(), probdesc::filename, first_invar(), gtopology::gelements, mechtop::give_ndofn(), mechtop::give_node_code_numbers(), lhsrhs::give_rhs(), mechmat::givestress(), Gtm, mechmat::ip, probdesc::loadcoef, Lsrs, Mm, Mp, Mt, intpoints::ncompstr, Ndofm, gtopology::ne, mechtop::ne, newton_raphsonep(), probdesc::niep, nod, gelement::nodes, nullv(), p, probdesc::path, print_close(), print_init(), springel::res_internal_forces(), Spring, intpoints::ssst, stress, mechmat::tnip, vector_tensor(), and write_elemscalar().
Referenced by epressure_solver().
| void general_epressure | ( | long | lcid | ) |
earth pressure solver with load iteration
The function solves erath pressure problem using GLPT. Iteration of load is used.
| lcid | - load case id. |
Created by Tomas Koudelka
Definition at line 96 of file epsolver.cpp.
References arclengthrv(), dloadcase::assemble(), loadcase::assemble(), mechbclc::dlc, probdesc::errep, f, lhsrhs::give_lhs(), lhsrhs::give_rhs(), mechbclc::lc, Lsrs, Mb, memset(), Mp, Ndofm, probdesc::niep, nullv(), print_close(), print_err(), print_init(), print_step(), probdesc::rlerrep, ss(), probdesc::stepep, and probdesc::zero.
Referenced by epressure_solver().
| void general_epressure_varsup | ( | long | lcid | ) |
earth pressure solver with temporary supports
The function solves erath pressure problem using GLPT. Temporary supports are used for better algorithm convergence.
| lcid | - load case id. |
Created by Tomas Koudelka
Definition at line 185 of file epsolver.cpp.
References lhsrhs::alloc(), arclengthrv(), dloadcase::assemble(), loadcase::assemble(), gnode::cn, mechbclc::dlc, probdesc::errep, f, gtopology::gencodnum(), lhsrhs::give_lhs(), gtopology::give_ndofn(), lhsrhs::give_rhs(), gtopology::gnodes, Gtm, internal_forces(), mechbclc::lc, Lsrs, Mb, memset(), Mp, Mt, Ndofm, probdesc::niep, gtopology::nn, mechtop::nn, noddispl(), probdesc::nselnodep, nullv(), print_close(), print_init(), print_step(), gtopology::restore_codnum(), restore_displ(), probdesc::rlerrep, probdesc::selnodep, Smat, ss(), and probdesc::zero.
Referenced by epressure_solver().
| double newton_raphsonep | ( | long | lcid, | |
| double * | pfp | |||
| ) |
the Newton-Raphsonmethop modified for starting from a reached state
This function solves system of nonlinear algebraic equations by the arc-length method. This method was modified to stop iteration on the given value rlambda of the lambda parameter. Rlambda value is reached with required error rlerr. Only one right hand side vector is supported with respect to nonlinearity and absence of superposition.
| lcid | - load case id | |
| rlambda | - required value of lambda parameter | |
| rlerr | - required error between lambda and rlambda |
double arclengthrv (long lcid, double rlambda, double rlerr)
d odpovida delte dd odpovida DELTE (trojuhelniku)
fc - konstantni vektor fp - vektor proporcionality
n - pocet radku matice { long i,j,k,n,ni,ini,stop,modif,li, kk; double a0,a1,a2,d,l1,l2,dl,dlmax,dlmin,psi,lambda,blambda,dlambda,ddlambda,norf,norfp,norfa,norv,zero,ierr; double lambdao,ss1,ss2,ss3,ss4,ss5; double *r,*ra,*ddr,*u,*v,*f,*fa,*fp,*fc,*fi; FILE *grout;
number of rows of the matrix n = Ndofm; maximum number of increments ni = Mp->nial; maximum number of iterations in one increment ini = Mp->niilal; computer zero zero = Mp->zero; required error in inner loop ierr = Mp->erral; length of arc dl = Mp->dlal; maximum length of arc dlmax = Mp->dlmaxal; minimum length of arc dlmin = Mp->dlminal; displacement-loading driving switch psi = Mp->psial;
allocation of auxiliary arrays r = new double [n]; ddr = new double [n]; u = new double [n]; v = new double [n]; f = new double [n]; fa = new double [n]; fi = new double [n]; fc = new double [n]; memset(fc, 0, sizeof(*fc)*n);
allocation of backup arrays for state variables in integration points Mm->alloc_backup ();
initialization phase ra = Lsrs->give_lhs (lcid); fp = Lsrs->give_rhs (lcid); fc = Lsrs->give_rhs (lcid*2+1);
array containing selected DOFs seldofinit ();
lambda=0.0; lambdao=0.0; li=0;
norm of proportionality vector norfp = ss (fp,fp,n); matrix sm(6,6); modif=0; grout = fopen("gr.dat", "wt");
main iteration loop ****
for (i=li;i<ni;i++){
fprintf(grout, "%e\n", lambda); for (kk = 0; kk < Mt->nn; kk++) fprintf(grout, " %e\n", ra[3*kk]); fprintf(grout, "Tuhosti --------\n"); for (kk = 15; kk < Mt->ne; kk++) { Spring->res_stiffness_matrix (kk,sm); fprintf(grout, " %e\n", sm[0][0]); } fprintf(grout, "\n"); fflush(grout); backup of left hand side vector for (j=0;j<n;j++) r[j]=ra[j]; backup of reached lambda parameter blambda=lambda; backup of state variables Mm->backup_state_var (); aux_nonlin_print (gr,ra,lambda);
fprintf (stdout,"\n arc-length: increment %ld lambda %e dl %e",i,lambda,dl); if (Mespr==1) { if (Mp->adp==1) Mp->ado->print (Out,i,lambda,0); }
assembling of tangent stiffness matrix stiffness_matrix (lcid);
backup of the fp, in ldl_sky the right hand side will be destroyed for (j=0;j<n;j++) f[j]=fp[j]; solution of K(r).v=F Smat->solve_system (v,f);
generalized norm of displacement increments norv = displincr (v,n);
dlambda = dl/sqrt(norv+psi*psi*norfp);
for (j=0;j<n;j++) { ddr[j]=dlambda*v[j]; ra[j]+=ddr[j]; fa[j]=fc[j]+(lambda+dlambda)*fp[j]; }
ddlambda=dlambda;
computation of internal forces internal_forces (lcid,fi);
for (k=0;k<n;k++) f[k]=fa[k]-fi[k];
norf=ss(f,f,n); norf=sqrt(norf); norfa = ss(fa,fa,n); norfa = sqrt(norfa); norf /= norfa;
if (Mespr==1) fprintf (stdout,"\n %e %e norf %e",lambda,dl,norf);
if (norf<ierr) {
no inner iteration loop is required ***
modif++; if (modif>1) { arc length modification dl*=2.0; if (dl>dlmax) dl=dlmax; modif=0;
if (Mespr==1) { fprintf (stdout,"\n arc length must be modified (dl increase) dl=%e norf=%e",dl,norf); fprintf (Out,"\n arc length must be modified (dl increase) dl=%e norf=%e",dl,norf); } }
lambda+=dlambda; if (((lambda - rlambda) > 0.0) && (fabs(lambda-rlambda) > rlerr)) { modification of the arc length dl/=2.0; if (dl<dlmin) { dl=dlmin; break; }
if (Mespr==1) { fprintf (stdout,"\n arc length must be modified (dl decrease) dl=%e norf=%e",dl,norf); fprintf (Out,"\n arc length must be modified (dl decrease) dl=%e norf=%e",dl,norf); }
restoring of left hand side vector
for (j=0;j<n;j++) ra[j]=r[j];
restoring of lambda parameter lambda=blambda; } if (fabs(lambda-rlambda) <= rlerr) break; Mm->updateipval(); continue; }
inner iteration loop ****
stop=0; for (j=0;j<ini;j++) { Smat->solve_system (u,f); Mp->ssle.solve_system (Gtm,Smat,u,f,Out);
for (k=0;k<n;k++) { f[k]=ddr[k]; ddr[k]+=u[k]; }
coefficient of quadratic equation a2=norv+psi*psi*norfp; a1=2.0*ss(ddr,v,n)+2.0*ddlambda*psi*psi*norfp; a0=ss(ddr,ddr,n)+ddlambda*ddlambda*psi*psi*norfp-dl*dl;
solution of quadratic equation if (fabs(a2)<zero) { if (fabs(a1)<zero) { if (fabs(a0)>zero) fprintf (stderr,"\n\n nonsolvable constrained condition in function arclength"); else fprintf (stderr,"\n\n infinite number of solution of constrained condition in function arclength"); } else dlambda=(0.0-a0)/a1; } else { d=a1*a1-4.0*a2*a0; if (d<0.0) { fprintf (stderr,"\n negative discriminant in function arclength"); break; } l1=(0.0-a1+sqrt(d))/2.0/a2; l2=(0.0-a1-sqrt(d))/2.0/a2; }
zacatek novinky
ss1=0.0; ss2=0.0; ss3=0.0; ss4=0.0; ss5=0.0; for (k=0;k<n;k++) { ss1+=(ddr[k]+l1*v[k])*f[k]; ss2+=(ddr[k]+l2*v[k])*f[k]; ss3+=(ddr[k]+l1*v[k])*(ddr[k]+l1*v[k]); ss4+=(ddr[k]+l2*v[k])*(ddr[k]+l2*v[k]); ss5+=f[k]*f[k]; }
if (ss1/ss3/ss5>ss2/ss4/ss5) dlambda=l1; else dlambda=l2; konec novinky fprintf(Out, "Dlambda %g, i = %ld, j = %ld\n", dlambda, i, j);
for (k=0;k<n;k++) { ddr[k]+=dlambda*v[k]; ra[k]+=u[k]+dlambda*v[k]; fa[k]+=dlambda*fp[k]; } ddlambda+=dlambda;
computation of internal forces internal_forces (lcid,fi);
for (k=0;k<n;k++) f[k]=fa[k]-fi[k];
norf=ss(f,f,n); norf=sqrt(norf); norfa = ss(fa,fa,n); norfa=sqrt(norfa); norf /= norfa;
if (Mespr==1) { fprintf (stdout,"\n increment %ld inner loop %ld norf=%e",i,j,norf); fprintf (Out,"\n increment %ld inner loop %ld norf=%e",i,j,norf); }
if (norf<ierr) { lambda+=ddlambda;
stop=1; if (((lambda - rlambda) > 0.0) && (fabs(lambda-rlambda) > rlerr)) stop=0; if (fabs(lambda-rlambda) <= rlerr) stop=2; if (stop > 0) Mm->updateipval(); break; } }
limit value of lambda was reached if (stop == 2) break;
modif=0;
if (stop==0) { modification of the arc length dl/=2.0; if (dl<dlmin) { dl=dlmin; break; }
if (Mespr==1) { fprintf (stdout,"\n arc length must be modified (dl decrease) dl=%e norf=%e",dl,norf); fprintf (Out,"\n arc length must be modified (dl decrease) dl=%e norf=%e",dl,norf); }
restoring of left hand side vector
for (j=0;j<n;j++) ra[j]=r[j];
restoring of lambda parameter lambda=blambda; restoring of state variables Mm->restore_state_var (); } } fclose(grout); delete [] r; delete [] fi; delete [] fa; delete [] fc; delete [] f; delete [] v; delete [] u; delete [] ddr; return lambda; } Function solves system of nonlinear algebraic equations by Newton-Raphson method. Solved system does not contain time variable, it takes into account previous reached displacement/stress state and starts with nonzero displacement vector. pfp is proportional load vector of the previous reached initial displacements.
| lcid | - loadcase id | |
| pfp | - proportional load vector of the prvious reached initial displacements. |
The function returns reached load coefficient lambda.
Created by Tomas Koudelka, 16.8.2001
Definition at line 1198 of file epsolver.cpp.
References nonlinman::errnr, f, lhsrhs::give_lhs(), lhsrhs::give_rhs(), Gtm, nonlinman::incrnr, internal_forces(), Lsrs, memset(), Mespr, nonlinman::minincrnr, Mm, Mp, Ndofm, nonlinman::niilnr, nonlinman::ninr, probdesc::nlman, Out, print_flush(), print_step(), Smat, slesolv::solve_system(), ss(), probdesc::ssle, stiffness_matrix(), mechmat::updateipval(), and probdesc::zero.
Referenced by femplast_epressure().
| void restore_displ | ( | long | lcid, | |
| double ** | bckr | |||
| ) |
extraction of the reached displacement form the backup vector
The function extracts nodal displacements bckr to the displacement vector of Lsrs.
| lcid | - number of load case | |
| bckr | - array with displacement in the each node |
Created by Tomas Koudelka 3.12.2002
Definition at line 436 of file epsolver.cpp.
References gnode::cn, gtopology::gnodes, Gtm, lhsrhs::lhs, Lsrs, Ndofm, gnode::ndofn, and gtopology::nn.
Referenced by general_epressure_varsup().
| void solve_epressure | ( | ) |
main computation function for earth pressure problem type
The function controls run of earth pressure problem computation.
| The | function does not return anything. |
Created by Tomas Koudelka
Definition at line 29 of file epsolver.cpp.
References dloadcase::assemble(), loadcase::assemble(), mechbclc::dlc, epplast_sol, epressure_solver(), lhsrhs::give_lhs(), lhsrhs::give_rhs(), mechbclc::lc, Lsrs, Mb, Mp, Ndofm, lhsrhs::nlc, nullv(), and probdesc::tepsol.
Referenced by solve_mefel_deterministic_problem().
1.6.2