Changeset 77

Timestamp:

04/28/09 12:08:24 (16 years ago)

Author:

Mathieu Morlighem

Message:

Fixed differences between ice and cielo (MOUNTAINEXPONENT applied on alpha2 and not K)

Location:

issm/trunk/src

Files:

• c/objects/Tria.cpp (modified) (1 diff)
• m/classes/@triaelem/CreateKMatrix.m (modified) (6 diffs)

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

@@ -336 +336 @@
 
         double MAXSLOPE=.06; // 6 %
-        double MOUNTAINKEXPONENT=10;
+        double MOUNTAINKEXPONENT=12;
 
 

issm/trunk/src/m/classes/@triaelem/CreateKMatrix.m

@@ -207 +207 @@
 %some parameters
 MAXSLOPE=.06;  %any element with slope>MAXSLOPE is considered a "rock" element, with infinite stiffness.
-MOUNTAINKEXPONENT=10; % "infinite" stiffness is going to be  10^MOUNTAINKEXPONENT
+MOUNTAINKEXPONENT=12; % "infinite" stiffness is going to be  10^MOUNTAINKEXPONENT
 
 %Create elementary stiffness matrix 
@@ -323 +323 @@
         frictionparameters.k=K_list;
 
-        alpha2=Getalpha2(frictionparameters);
+        alpha2_list=Getalpha2(frictionparameters);
         if (element_debug & triaelem.id==element_debugid),
                 disp(sprintf('   alpha2_list (%g %g %g )',alpha2_list(1),alpha2_list(2),alpha2_list(3)));
@@ -368 +368 @@
 
                 if (slope_magnitude>MAXSLOPE),
-                        K=10^(-MOUNTAINKEXPONENT);
+                        alpha2_list=10^(MOUNTAINKEXPONENT)*ones(3,1);
                 end
         end
@@ -398 +398 @@
                 matice.B=B_param; clear B_param.
         end
-
 
         %Get viscosity at last two iterations: 
@@ -448 +447 @@
  
                 %compute alpha2 for the current gaussian point
-                alpha2_g=GetParameterValue(triaelem,alpha2,gauss_l1l2l3);
+                alpha2_g=GetParameterValue(triaelem,alpha2_list,gauss_l1l2l3);
         
                 if velocity_is_present
@@ -474 +473 @@
                 L
         end
-
 end %for ig=1:num_gauss,
-        
+
 if (element_debug & triaelem.id==element_debugid),
         disp(sprintf('Ke_gg->terms:\n'));