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

../trunk-jpl/src/c/solvers/solver_newton.cpp

@@ -55 +55 @@
                 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 +101 @@
                         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++;
         }