source: issm/oecreview/Archive/14064-14311/ISSM-14166-14167.diff@ 14312

Index: ../trunk-jpl/src/c/solvers/solver_newton.cpp
===================================================================
--- ../trunk-jpl/src/c/solvers/solver_newton.cpp        (revision 14166)
+++ ../trunk-jpl/src/c/solvers/solver_newton.cpp        (revision 14167)
@@ -55,12 +55,32 @@
                xdelete(&old_uf);old_uf=uf;
 
                /*Solver forward model*/
-               femmodel->SystemMatricesx(&Kff, &Kfs, &pf, &df, NULL);
+               if(count==1){
+                       femmodel->SystemMatricesx(&Kff, &Kfs, &pf, &df, NULL);
+                       CreateNodalConstraintsx(&ys,femmodel->nodes,configuration_type);
+                       Reduceloadx(pf,Kfs,ys);xdelete(&Kfs);
+                       Solverx(&uf,Kff,pf,old_uf,df,femmodel->parameters);xdelete(&df);
+                       Mergesolutionfromftogx(&ug,uf,ys,femmodel->nodes,femmodel->parameters);xdelete(&ys);
+                       InputUpdateFromSolutionx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,ug);xdelete(&ug);
+                       xdelete(&old_ug);old_ug=ug;
+                       xdelete(&old_uf);old_uf=uf;
+               }
+               uf=old_uf->Duplicate(); old_uf->Copy(uf);
+
+               /*Prepare next iteration using Newton's method*/
+               femmodel->SystemMatricesx(&Kff, &Kfs, &pf, &df, NULL);xdelete(&df);
                CreateNodalConstraintsx(&ys,femmodel->nodes,configuration_type);
-               Reduceloadx(pf,Kfs,ys);xdelete(&Kfs);
-               Solverx(&uf,Kff,pf,old_uf,df,femmodel->parameters);xdelete(&df);
+               Reduceloadx(pf,Kfs,ys);   xdelete(&Kfs);
+
+               pJf=pf->Duplicate();
+               Kff->MatMult(uf,pJf);// xdelete(&Kff);
+               pJf->Scale(-1.0); pJf->AXPY(pf,+1.0);     //xdelete(&pf);
+
+               femmodel->CreateJacobianMatrixx(&Jff,kmax);
+               Solverx(&duf,Jff,pJf,NULL,NULL,femmodel->parameters); xdelete(&Jff); xdelete(&pJf);
+               uf->AXPY(duf, 1.0); xdelete(&duf);
                Mergesolutionfromftogx(&ug,uf,ys,femmodel->nodes,femmodel->parameters);xdelete(&ys);
-               InputUpdateFromSolutionx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,ug);xdelete(&ug);
+               InputUpdateFromSolutionx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,ug);
 
                /*Check convergence*/
                convergence(&converged,Kff,pf,uf,old_uf,femmodel->parameters); 
@@ -81,20 +101,6 @@
                        break;
                }
 
-               /*Prepare next iteration using Newton's method*/
-               femmodel->SystemMatricesx(&Kff, &Kfs, &pf, &df, NULL);
-               CreateNodalConstraintsx(&ys,femmodel->nodes,configuration_type);
-               Reduceloadx(pf,Kfs,ys);   xdelete(&Kfs);
-
-               pJf=pf->Duplicate(); Kff->MatMult(uf,pJf); xdelete(&Kff);
-               pJf->Scale(-1.0); pJf->AXPY(pf,+1.0);     xdelete(&pf);
-
-               femmodel->CreateJacobianMatrixx(&Jff,kmax);
-               Solverx(&duf,Jff,pJf,NULL,NULL,femmodel->parameters); xdelete(&Jff); xdelete(&pJf);
-               uf->AXPY(duf, 1.0); xdelete(&duf);
-               Mergesolutionfromftogx(&ug,uf,ys,femmodel->nodes,femmodel->parameters);xdelete(&ys);
-               InputUpdateFromSolutionx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,ug);
-
                count++;
        }