00001 #include <math.h>
00002 #include <string.h>
00003 #include <stdlib.h>
00004 #include "pglobal.h"
00005 #include "pnewtonraph.h"
00006 #include "nssolver.h"
00007 #include "backupsol.h"
00008 #include "global.h"
00009 #include "globmat.h"
00010 #include "mechprint.h"
00011 #include "gmatrix.h"
00012 #include "gtopology.h"
00013 #include "mechprint.h"
00014 #include "mathem.h"
00015 #include "vector.h"
00016 #include "matrix.h"
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033
00034
00035
00036
00037
00038
00039
00040
00041 double par_gnewton_raphson_one_step(long lcid, nonlinman *nlman, double *fa, double *ra, double *fb, double *dr, double *fi,
00042 long istep, long &j, long li, long ini, double ierr)
00043 {
00044 long n = Ndofm;
00045 double norf, norfa;
00046 matrix lsm_a(3,3);
00047 vector lsm_r(3), lsm_l(3);
00048
00049
00050 assemble_stiffness_matrix(lcid,istep,-1,li);
00051
00052
00053 Psolm->par_linear_solver (Gtm,Smat,dr,fb,Out,Mespr);
00054
00055
00056
00057 addv(ra, dr, n);
00058
00059
00060 internal_forces (lcid,fi);
00061
00062
00063
00064 subv(fa, fi, fb, n);
00065
00066
00067 norfa=Psolm->pss (fa,fa,Out);
00068
00069 norf=Psolm->pss(fb,fb,Out);
00070 if (norfa != 0.0)
00071 norf /= norfa;
00072
00073 if ((Myrank==0)&&(Mespr==1))
00074 fprintf (stdout,", norm %e",norf);
00075
00076 if (norf<ierr)
00077 {
00078 j = -1;
00079 return norf;
00080 }
00081
00082
00083 fillm(0.0, lsm_a);
00084 fillv(0.0, lsm_r);
00085 for (j=0; j<ini; j++)
00086 {
00087 Mp->jstep = j;
00088
00089
00090 assemble_stiffness_matrix(lcid,istep,j,li);
00091
00092
00093 Psolm->par_linear_solver (Gtm,Smat,dr,fb,Out,Mespr);
00094
00095
00096 addv(ra, dr, n);
00097
00098
00099 internal_forces(lcid,fi);
00100
00101
00102
00103 subv(fa, fi, fb, n);
00104
00105
00106 norf=Psolm->pss(fb,fb,Out);
00107 if (norfa != 0.0)
00108 norf /= norfa;
00109
00110 if ((Myrank==0) && (Mespr==1))
00111 fprintf (stdout,"\n increment %ld inner loop j %ld norf %e",istep,j,norf);
00112
00113 if (norf < ierr)
00114
00115 return norf;
00116
00117
00118 if (check_divergency(nlman, lsm_a, lsm_r, lsm_l, j, norf))
00119 return norf;
00120 }
00121
00122
00123 return norf;
00124 }