Changeset 24070

Timestamp:

07/09/19 05:22:14 (6 years ago)

Author:

Mathieu Morlighem

Message:

BUG: bad calculation of tx and ty

Location:

issm/trunk-jpl/src/c/classes/Loads

Files:

• Channel.cpp (modified) (6 diffs)
• Channel.h (modified) (2 diffs)

issm/trunk-jpl/src/c/classes/Loads/Channel.cpp

@@ -48 +48 @@
 
         /*Set channel cross section to 0*/
-        this->S = 0.;
+        this->S    = 0.;
+        this->Sold = 0.;
 
         /*Get edge info*/
@@ -394 +395 @@
         IssmDouble tx = xyz_list_tria[index2][0] - xyz_list_tria[index1][0];
         IssmDouble ty = xyz_list_tria[index2][1] - xyz_list_tria[index1][1];
-        tx = tx/sqrt(tx*tx+ty*ty);
-        ty = ty/sqrt(tx*tx+ty*ty);
+        IssmDouble Lt = sqrt(tx*tx+ty*ty);
+        tx = tx/Lt;
+        ty = ty/Lt;
 
         /* Start  looping on the number of gaussian points: */
@@ -524 +526 @@
         IssmDouble tx = xyz_list_tria[index2][0] - xyz_list_tria[index1][0];
         IssmDouble ty = xyz_list_tria[index2][1] - xyz_list_tria[index1][1];
-        tx = tx/sqrt(tx*tx+ty*ty);
-        ty = ty/sqrt(tx*tx+ty*ty);
+        IssmDouble Lt = sqrt(tx*tx+ty*ty);
+        tx = tx/Lt;
+        ty = ty/Lt;
 
         /* Start  looping on the number of gaussian points: */
@@ -588 +591 @@
 }
 /*}}}*/
+void           Channel::SetChannelCrossSectionOld(void){/*{{{*/
+
+        this->Sold = this->S;
+
+} /*}}}*/
 void           Channel::UpdateChannelCrossSection(void){/*{{{*/
 
@@ -630 +638 @@
         IssmDouble tx = xyz_list_tria[index2][0] - xyz_list_tria[index1][0];
         IssmDouble ty = xyz_list_tria[index2][1] - xyz_list_tria[index1][1];
-        tx = tx/sqrt(tx*tx+ty*ty);
-        ty = ty/sqrt(tx*tx+ty*ty);
+        IssmDouble Lt = sqrt(tx*tx+ty*ty);
+        tx = tx/Lt;
+        ty = ty/Lt;
 
         /*Evaluate fields on center of edge*/
@@ -681 +690 @@
 
         /*Solve ODE*/
-        this->S = ODE1(alpha,beta,this->S,dt,2);
+        this->S = ODE1(alpha,beta,this->Sold,dt,1);
 
         /*Make sure Area > 0*/

issm/trunk-jpl/src/c/classes/Loads/Channel.h

@@ -20 +20 @@
         private: 
                 IssmDouble S;
+                IssmDouble Sold;
                 bool       boundary;
 
@@ -76 +77 @@
                 /*}}}*/
                 /*Channel management:{{{*/
+                void           SetChannelCrossSectionOld(void);
                 void           UpdateChannelCrossSection(void);
                 ElementVector* CreatePVectorHydrologyGlaDS(void);