Changeset 11337

Timestamp:

02/06/12 16:41:32 (13 years ago)

Author:

Mathieu Morlighem

Message:

Fixed Newton's method: do not change ys

Location:

issm/trunk-jpl/src

Files:

• c/solvers/solver_newton.cpp (modified) (1 diff)
• m/solvers/solver_newton.m (modified) (1 diff)

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

@@ -68 +68 @@
 
                 /*Prepare next iteration using Newton's method*/
-                SystemMatricesx(&Kff,NULL,&pf,NULL,&kmax,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters);
-                CreateJacobianMatrixx(&Jff,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,kmax);
-                CreateNodalConstraintsx(&ys,femmodel->nodes,configuration_type);VecScale(ys,0.);
+                SystemMatricesx(&Kff,&Kfs,&pf,NULL,&kmax,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters);
+                CreateNodalConstraintsx(&ys,femmodel->nodes,configuration_type);
+                Reduceloadx(pf,Kfs,ys);   MatFree(&Kfs);
+
                 VecDuplicate(pf,&pJf);
                 MatMultPatch(Kff,uf,pJf); MatFree(&Kff);
                 VecScale(pJf,-1.);
                 VecAXPY(pJf,+1.,pf);      VecFree(&pf);
-                Solverx(&duf,Jff,pJf,NULL,NULL,femmodel->parameters);
-                MatFree(&Jff);            VecFree(&pJf);
+
+                CreateJacobianMatrixx(&Jff,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,kmax);
+                Solverx(&duf,Jff,pJf,NULL,NULL,femmodel->parameters); MatFree(&Jff);//VecFree(&pJf);
                 VecAXPY(uf,1.,duf);      VecFree(&duf);
                 Mergesolutionfromftogx(&ug,uf,ys,femmodel->nodes,femmodel->parameters);VecFree(&ys);

issm/trunk-jpl/src/m/solvers/solver_newton.m

@@ -51 +51 @@
                 Jff=CreateJacobianMatrix(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters,kmax);
                 duf=Solver(Jff,pJf,[],[],femmodel.parameters);
-                ys=0*ys;
                 uf=uf+duf;
                 ug=Mergesolutionfromftog(uf,ys,femmodel.nodes,femmodel.parameters);