Changeset 926

Timestamp:

06/11/09 16:45:42 (16 years ago)

Author:

Mathieu Morlighem

Message:

fixed melting

File:

• issm/trunk/src/c/objects/Pengrid.cpp (modified) (6 diffs)

issm/trunk/src/c/objects/Pengrid.cpp

@@ -309 +309 @@
 void  Pengrid::PenaltyCreateKMatrixMelting(Mat Kgg,void* vinputs,double kmax,int analysis_type,int sub_analysis_type){
 
+
         int found=0;
         const int numgrids=1;
@@ -325 +326 @@
         ParameterInputs* inputs=NULL;
 
+        /*check that pengrid is not a clone (penalty to be added only once)*/
+        if (node->IsClone()) return;
+
         /*recover pointers: */
         inputs=(ParameterInputs*)vinputs;
@@ -345 +349 @@
         if (temperature<t_pmp){ //If T<Tpmp, there must be no melting. Therefore, melting should be  constrained to 0 when T<Tpmp, instead of using spcs, use penalties
                 Ke[0][0]=kmax*pow(10,penalty_offset);
+                if (id==65){
+                        printf("Ke_terms[0] = %10.10lf\n",Ke[0]);
+                        printf("penalty_offset(K) = %g\n",penalty_offset);
+                        printf("kmax = %10.10lf\n",kmax);
+                }
         }
         
@@ -458 +467 @@
         ParameterInputs* inputs=NULL;
 
+        /*check that pengrid is not a clone (penalty to be added only once)*/
+        if (node->IsClone()) return;
+
         /*recover pointers: */
         inputs=(ParameterInputs*)vinputs;
@@ -495 +507 @@
                 if (sub_analysis_type==SteadyAnalysisEnum()){
                         P_terms[0]=melting_offset*pow(10,penalty_offset)*(temperature-t_pmp);
+                        if (id==65){
+                                printf("penalty_offset(P) = %g\n",penalty_offset);
+                                printf("melting_offset(P) = %10.10lf\n",melting_offset);
+                                printf("P_terms[0] = %10.10lf\n",P_terms[0]);
+                        }
                 }
                 else{
@@ -553 +570 @@
         ParameterInputs* inputs=NULL;
 
+        /*check that pengrid is not a clone (penalty to be added only once)*/
+        if (node->IsClone()){
+                unstable=0;
+                *punstable=unstable;
+                return;
+        }
+
         /*recover pointers: */
         inputs=(ParameterInputs*)vinputs;
-
 
         //First recover beta, pressure and temperature vectors;