Changeset 14415

Timestamp:

02/28/13 08:59:05 (12 years ago)

Author:

Mathieu Morlighem

Message:

NEW: reimplemented newton 2, where deltav is applied to the stiffness matrix

File:

• issm/trunk-jpl/src/c/solvers/solver_newton.cpp (modified) (4 diffs)

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

@@ -15 +15 @@
         /*intermediary: */
         bool   converged;
-        int    count,newton;
+        int    count;
         IssmDouble kmax;
         Matrix<IssmDouble>* Kff = NULL;
@@ -37 +37 @@
         femmodel->parameters->FindParam(&max_nonlinear_iterations,DiagnosticMaxiterEnum);
         femmodel->parameters->FindParam(&configuration_type,ConfigurationTypeEnum);
-        femmodel->parameters->FindParam(&newton,DiagnosticIsnewtonEnum);
         femmodel->UpdateConstraintsx();
 
@@ -57 +56 @@
 
                 /*Solver forward model*/
-                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);
+                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);
+                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);
+                InputUpdateFromSolutionx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,ug);xdelete(&ug);
 
                 /*Check convergence*/
@@ -103 +82 @@
                 }
 
+                /*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++;
         }