Changeset 11327

Timestamp:

02/06/12 08:50:20 (13 years ago)

Author:

Mathieu Morlighem

Message:

Added penalties in Newton's method

Location:

issm/trunk-jpl/src/c

Files:

• modules/CreateJacobianMatrixx/CreateJacobianMatrixx.cpp (modified) (2 diffs)
• modules/CreateJacobianMatrixx/CreateJacobianMatrixx.h (modified) (1 diff)
• objects/Loads/Icefront.h (modified) (1 diff)
• objects/Loads/Load.h (modified) (1 diff)
• objects/Loads/Numericalflux.h (modified) (1 diff)
• objects/Loads/Pengrid.h (modified) (1 diff)
• objects/Loads/Penpair.cpp (modified) (2 diffs)
• objects/Loads/Penpair.h (modified) (1 diff)
• objects/Loads/Riftfront.h (modified) (1 diff)
• objects/Numerics/ElementMatrix.cpp (modified) (2 diffs)
• solvers/solver_newton.cpp (modified) (4 diffs)

issm/trunk-jpl/src/c/modules/CreateJacobianMatrixx/CreateJacobianMatrixx.cpp

@@ -10 +10 @@
 #include "../../EnumDefinitions/EnumDefinitions.h"
 
-void CreateJacobianMatrixx(Mat* pJff,Elements* elements,Nodes* nodes, Vertices* vertices,Loads* loads,Materials* materials, Parameters* parameters){
+void CreateJacobianMatrixx(Mat* pJff,Elements* elements,Nodes* nodes, Vertices* vertices,Loads* loads,Materials* materials, Parameters* parameters,double kmax){
         
-        int      connectivity;
+        int      i,connectivity;
         int      numberofdofspernode;
         int      fsize,configuration_type;
         Element *element = NULL;
+        Load    *load    = NULL;
         Mat      Jff     = NULL;
 
@@ -31 +32 @@
         
         /*Create and assemble matrix*/
-        for(int i=0;i<elements->Size();i++){
+        for(i=0;i<elements->Size();i++){
                 element=(Element*)elements->GetObjectByOffset(i);
                 element->CreateJacobianMatrix(Jff);
+        }
+        for (i=0;i<loads->Size();i++){
+                load=(Load*)loads->GetObjectByOffset(i);
+                if(load->InAnalysis(configuration_type)) load->CreateJacobianMatrix(Jff);
+                if(load->InAnalysis(configuration_type)) load->PenaltyCreateJacobianMatrix(Jff,kmax);
         }
         MatAssemblyBegin(Jff,MAT_FINAL_ASSEMBLY);

issm/trunk-jpl/src/c/modules/CreateJacobianMatrixx/CreateJacobianMatrixx.h

@@ -10 +10 @@
 
 /* local prototypes: */
-void CreateJacobianMatrixx(Mat* pJff,Elements* elements,Nodes* nodes, Vertices* vertices,Loads* loads,Materials* materials, Parameters* parameters);
+void CreateJacobianMatrixx(Mat* pJff,Elements* elements,Nodes* nodes, Vertices* vertices,Loads* loads,Materials* materials, Parameters* parameters,double kmax);
 
 #endif  /* _CREATEJACOBIANMATRIXX_H */

issm/trunk-jpl/src/c/objects/Loads/Icefront.h

@@ -76 +76 @@
                 void  CreateKMatrix(Mat Kff, Mat Kfs);
                 void  CreatePVector(Vec pf);
+                void  CreateJacobianMatrix(Mat Jff){_error_("Not implemented yet");};
+                void  PenaltyCreateJacobianMatrix(Mat Jff,double kmax){_error_("Not implemented yet");};
                 void  PenaltyCreateKMatrix(Mat Kff, Mat kfs, double kmax);
                 void  PenaltyCreatePVector(Vec pf, double kmax);

issm/trunk-jpl/src/c/objects/Loads/Load.h

@@ -29 +29 @@
                 virtual void  CreateKMatrix(Mat Kff, Mat Kfs)=0;
                 virtual void  CreatePVector(Vec pf)=0;
+                virtual void  CreateJacobianMatrix(Mat Jff)=0;
+                virtual void  PenaltyCreateJacobianMatrix(Mat Jff,double kmax)=0;
                 virtual void  PenaltyCreateKMatrix(Mat Kff, Mat Kfs, double kmax)=0;
                 virtual void  PenaltyCreatePVector(Vec pf, double kmax)=0;

issm/trunk-jpl/src/c/objects/Loads/Numericalflux.h

@@ -72 +72 @@
                 void  CreateKMatrix(Mat Kff, Mat Kfs);
                 void  CreatePVector(Vec pf);
+                void  CreateJacobianMatrix(Mat Jff){_error_("Not implemented yet");};
+                void  PenaltyCreateJacobianMatrix(Mat Jff,double kmax){_error_("Not implemented yet");};
                 void  PenaltyCreateKMatrix(Mat Kff, Mat kfs, double kmax);
                 void  PenaltyCreatePVector(Vec pf, double kmax);

issm/trunk-jpl/src/c/objects/Loads/Pengrid.h

@@ -77 +77 @@
                 void  CreateKMatrix(Mat Kff, Mat Kfs);
                 void  CreatePVector(Vec pf);
+                void  CreateJacobianMatrix(Mat Jff){_error_("Not implemented yet");};
+                void  PenaltyCreateJacobianMatrix(Mat Jff,double kmax){_error_("Not implemented yet");};
                 void  PenaltyCreateKMatrix(Mat Kff, Mat kfs, double kmax);
                 void  PenaltyCreatePVector(Vec pf, double kmax);

issm/trunk-jpl/src/c/objects/Loads/Penpair.cpp

@@ -209 +209 @@
 void  Penpair::CreatePVector(Vec pf){
 
+        /*No loads applied, do nothing: */
+        return;
+
+}
+/*}}}1*/
+/*FUNCTION Penpair::CreateJacobianMatrix{{{1*/
+void  Penpair::CreateJacobianMatrix(Mat Jff){
+        /*If you code this piece, don't forget that a penalty will be inactive if it is dealing with clone nodes*/
         /*No loads applied, do nothing: */
         return;
@@ -244 +252 @@
         /*No loads applied, do nothing: */
         return;
+}
+/*}}}1*/
+/*FUNCTION Penpair::PenaltyCreateJacobianMatrix{{{1*/
+void  Penpair::PenaltyCreateJacobianMatrix(Mat Jff,double kmax){
+        this->PenaltyCreateKMatrix(Jff,NULL,kmax);
 }
 /*}}}1*/

issm/trunk-jpl/src/c/objects/Loads/Penpair.h

@@ -64 +64 @@
                 void  CreateKMatrix(Mat Kff, Mat Kfs);
                 void  CreatePVector(Vec pf);
-                void  PenaltyCreateKMatrix(Mat Kff, Mat kfs, double kmax);
+                void  CreateJacobianMatrix(Mat Jff);
+                void  PenaltyCreateKMatrix(Mat Kff,Mat Kfs,double kmax);
                 void  PenaltyCreatePVector(Vec pf, double kmax);
+                void  PenaltyCreateJacobianMatrix(Mat Jff,double kmax);
                 bool  InAnalysis(int analysis_type);
                 /*}}}*/

issm/trunk-jpl/src/c/objects/Loads/Riftfront.h

@@ -84 +84 @@
                 void  CreateKMatrix(Mat Kff, Mat Kfs);
                 void  CreatePVector(Vec pf);
+                void  CreateJacobianMatrix(Mat Jff){_error_("Not implemented yet");};
+                void  PenaltyCreateJacobianMatrix(Mat Jff,double kmax){_error_("Not implemented yet");};
                 void  PenaltyCreateKMatrix(Mat Kff, Mat kfs, double kmax);
                 void  PenaltyCreatePVector(Vec pf, double kmax);

issm/trunk-jpl/src/c/objects/Numerics/ElementMatrix.cpp

@@ -251 +251 @@
         double* localvalues=NULL;
 
+        /*If Kfs is not provided, call the other function*/
+        if(!Kfs){
+                this->AddToGlobal(Kff);
+                return;
+        }
+
         /*In debugging mode, check consistency (no NaN, and values not too big)*/
         this->CheckConsistency();
@@ -300 +306 @@
         double* localvalues=NULL;
 
+        /*Check that Jff is not NULL*/
+        _assert_(Jff); 
+
         /*In debugging mode, check consistency (no NaN, and values not too big)*/
         this->CheckConsistency();

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

@@ -14 +14 @@
 
         /*intermediary: */
+        bool   converged;
+        int    num_unstable_constraints;
+        int    count;
+        double kmax;
         Mat Kff = NULL, Kfs    = NULL, Jff = NULL;
         Vec ug  = NULL, old_ug = NULL;
@@ -20 +24 @@
         Vec df  = NULL;
         Vec ys  = NULL;
-        
-        bool converged;
-        int num_unstable_constraints;
-        int count;
 
         /*parameters:*/
@@ -69 +69 @@
 
                 /*Prepare next iteration using Newton's method*/
-                CreateJacobianMatrixx(&Jff,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters);
-                SystemMatricesx(&Kff,NULL,&pf,NULL,NULL,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters);
+                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);
                 VecDuplicate(pf,&pJf);
@@ -78 +78 @@
                 Solverx(&duf,Jff,pJf,NULL,NULL,femmodel->parameters);
                 MatFree(&Jff);            VecFree(&pJf);
-                VecAXPY(uf,+1.,duf);      VecFree(&duf);
+                VecAXPY(uf,1.,duf);      VecFree(&duf);
                 Mergesolutionfromftogx(&ug,uf,ys,femmodel->nodes,femmodel->parameters);VecFree(&ys);
                 InputUpdateFromSolutionx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,ug);