Changeset 26294

Timestamp:

05/28/21 14:44:26 (4 years ago)

Author:

tsantos

Message:

CHG: updated viscosity computation for MLHO: now we have 4 expressions obtained by numerical integration (it is more consistent).

Location:

issm/trunk-jpl/src/c

Files:

• analyses/StressbalanceAnalysis.cpp (modified) (3 diffs)
• classes/Elements/Element.cpp (modified) (2 diffs)
• classes/Elements/Element.h (modified) (1 diff)
• classes/Materials/Material.h (modified) (1 diff)
• classes/Materials/Matestar.cpp (modified) (1 diff)
• classes/Materials/Matestar.h (modified) (1 diff)
• classes/Materials/Matice.cpp (modified) (1 diff)
• classes/Materials/Matice.h (modified) (1 diff)
• classes/Materials/Matlitho.h (modified) (1 diff)

issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp ¶

@@ -2882 +2882 @@
 
         /*Intermediaries*/
-        IssmDouble  viscosity,Jdet,thickness,n;
+        IssmDouble  viscosity[4],Jdet,thickness,n;
         IssmDouble *xyz_list = NULL;
         int      domaintype;
@@ -2917 +2917 @@
         Input* vx_input       = element->GetInput(VxEnum);        _assert_(vx_input); //vertically integrated vx
         Input* vy_input       = element->GetInput(VyEnum);        _assert_(vy_input); //vertically integrated vy
-        Input* vxshear_input  = element->GetInput(VxShearEnum);   _assert_(vxshear_input); 
-        Input* vyshear_input  = element->GetInput(VyShearEnum);   _assert_(vyshear_input);
+        Input* vxbase_input   = element->GetInput(VxBaseEnum);    _assert_(vxbase_input); 
+        Input* vybase_input   = element->GetInput(VyBaseEnum);    _assert_(vybase_input);
+        Input* vxshear_input  = element->GetInput(VxShearEnum);   _assert_(vxshear_input); //shear vx at the surface 
+        Input* vyshear_input  = element->GetInput(VyShearEnum);   _assert_(vyshear_input); //shear vy at the surface
         Input* n_input                   = element->GetInput(MaterialsRheologyNEnum); _assert_(n_input);
 
@@ -2934 +2936 @@
                 n_input->GetInputValue(&n,gauss);
                 int dim=2;//itapopo
-                element->material->ViscosityMLHO(&viscosity,dim,xyz_list,gauss,vx_input,vy_input,vxshear_input,vyshear_input,thickness_input);
-                //viscosity=10e13;//itapopo
+                element->material->ViscosityMLHO(&viscosity[0],dim,xyz_list,gauss,vxbase_input,vybase_input,vxshear_input,vyshear_input,thickness_input,n_input);
 
                 for(int i=0;i<numnodes;i++){//shape functions on tria element
                         for(int j=0;j<numnodes;j++){
-                                Ke->values[(4*i+0)*2*2*numnodes+4*j+0] += gauss->weight*Jdet*viscosity*thickness*(
+                                Ke->values[(4*i+0)*2*2*numnodes+4*j+0] += gauss->weight*Jdet*viscosity[0]*thickness*(
                                                         4.*dbasis[0*numnodes+j]*dbasis[0*numnodes+i] + dbasis[1*numnodes+j]*dbasis[1*numnodes+i]
                                                         );//K11
-                                Ke->values[(4*i+0)*2*2*numnodes+4*j+1] += gauss->weight*Jdet*viscosity*thickness*(
+                                Ke->values[(4*i+0)*2*2*numnodes+4*j+1] += gauss->weight*Jdet*viscosity[1]*thickness*(
                                                         4.*dbasis[0*numnodes+j]*dbasis[0*numnodes+i] + dbasis[1*numnodes+j]*dbasis[1*numnodes+i]
                                                         )*(n+1)/(n+2);//K12
-                                Ke->values[(4*i+0)*2*2*numnodes+4*j+2] += gauss->weight*Jdet*viscosity*thickness*(
+                                Ke->values[(4*i+0)*2*2*numnodes+4*j+2] += gauss->weight*Jdet*viscosity[0]*thickness*(
                                                         2.*dbasis[1*numnodes+j]*dbasis[0*numnodes+i] + dbasis[0*numnodes+j]*dbasis[1*numnodes+i]
                                                         );//K13
-                                Ke->values[ (4*i+0)*2*2*numnodes+4*j+3] += gauss->weight*Jdet*viscosity*thickness*(
+                                Ke->values[ (4*i+0)*2*2*numnodes+4*j+3] += gauss->weight*Jdet*viscosity[1]*thickness*(
                      2.*dbasis[1*numnodes+j]*dbasis[0*numnodes+i] + dbasis[0*numnodes+j]*dbasis[1*numnodes+i]
                      )*(n+1)/(n+2);//K14
                                 
-                                Ke->values[(4*i+1)*2*2*numnodes+4*j+0] += gauss->weight*Jdet*viscosity*thickness*(
+                                Ke->values[(4*i+1)*2*2*numnodes+4*j+0] += gauss->weight*Jdet*viscosity[1]*thickness*(
                                                         4.*dbasis[0*numnodes+j]*dbasis[0*numnodes+i] + dbasis[1*numnodes+j]*dbasis[1*numnodes+i]
                                                         )*(n+1)/(n+2);//K21
-                                Ke->values[(4*i+1)*2*2*numnodes+4*j+1] += gauss->weight*Jdet*viscosity*thickness*(
+                                Ke->values[(4*i+1)*2*2*numnodes+4*j+1] += gauss->weight*Jdet*viscosity[2]*thickness*(
                                                         4.*dbasis[0*numnodes+j]*dbasis[0*numnodes+i] + dbasis[1*numnodes+j]*dbasis[1*numnodes+i]
                                                         )*2*pow(n+1,2)/((2*n+3)*(n+2)) 
                                                         + 
-                                                        gauss->weight*Jdet*viscosity*basis[j]*basis[i]*pow(n+1,2)/(thickness*(2*n+1))
+                                                        gauss->weight*Jdet*viscosity[3]*basis[j]*basis[i]*pow(n+1,2)/(thickness*(2*n+1))
                                                         ;//K22
-                                Ke->values[(4*i+1)*2*2*numnodes+4*j+2] += gauss->weight*Jdet*viscosity*thickness*(
+                                Ke->values[(4*i+1)*2*2*numnodes+4*j+2] += gauss->weight*Jdet*viscosity[1]*thickness*(
                                                         2.*dbasis[1*numnodes+j]*dbasis[0*numnodes+i] + dbasis[0*numnodes+j]*dbasis[1*numnodes+i]
                                                         )*(n+1)/(n+2);//K23
-                                Ke->values[(4*i+1)*2*2*numnodes+4*j+3] += gauss->weight*Jdet*viscosity*thickness*(
+                                Ke->values[(4*i+1)*2*2*numnodes+4*j+3] += gauss->weight*Jdet*viscosity[2]*thickness*(
                                                         2.*dbasis[1*numnodes+j]*dbasis[0*numnodes+i] + dbasis[0*numnodes+j]*dbasis[1*numnodes+i]
                                                         )*2*pow(n+1,2)/((2*n+3)*(n+2));//K24
                                 
-                                Ke->values[(4*i+2)*2*2*numnodes+4*j+0] += gauss->weight*Jdet*viscosity*thickness*(
+                                Ke->values[(4*i+2)*2*2*numnodes+4*j+0] += gauss->weight*Jdet*viscosity[0]*thickness*(
                                                         2.*dbasis[0*numnodes+j]*dbasis[1*numnodes+i] + dbasis[1*numnodes+j]*dbasis[0*numnodes+i]
                                                         );//K31
-                                Ke->values[(4*i+2)*2*2*numnodes+4*j+1] += gauss->weight*Jdet*viscosity*thickness*(
+                                Ke->values[(4*i+2)*2*2*numnodes+4*j+1] += gauss->weight*Jdet*viscosity[1]*thickness*(
                                                         2.*dbasis[0*numnodes+j]*dbasis[1*numnodes+i] + dbasis[1*numnodes+j]*dbasis[0*numnodes+i]
                                                         )*(n+1)/(n+2);//K32
-                                Ke->values[(4*i+2)*2*2*numnodes+4*j+2] += gauss->weight*Jdet*viscosity*thickness*(
+                                Ke->values[(4*i+2)*2*2*numnodes+4*j+2] += gauss->weight*Jdet*viscosity[0]*thickness*(
                                                         4.*dbasis[1*numnodes+j]*dbasis[1*numnodes+i] + dbasis[0*numnodes+j]*dbasis[0*numnodes+i]
                                                         );//K33
-                                Ke->values[(4*i+2)*2*2*numnodes+4*j+3] += gauss->weight*Jdet*viscosity*thickness*(
+                                Ke->values[(4*i+2)*2*2*numnodes+4*j+3] += gauss->weight*Jdet*viscosity[1]*thickness*(
                                                         4.*dbasis[1*numnodes+j]*dbasis[1*numnodes+i] + dbasis[0*numnodes+j]*dbasis[0*numnodes+i]
                                                         )*(n+1)/(n+2);//K34
 
-                                Ke->values[(4*i+3)*2*2*numnodes+4*j+0] += gauss->weight*Jdet*viscosity*thickness*(
+                                Ke->values[(4*i+3)*2*2*numnodes+4*j+0] += gauss->weight*Jdet*viscosity[1]*thickness*(
                      2.*dbasis[0*numnodes+j]*dbasis[1*numnodes+i] + dbasis[1*numnodes+j]*dbasis[0*numnodes+i]
                      )*(n+1)/(n+2);//K41
-                                Ke->values[(4*i+3)*2*2*numnodes+4*j+1] += gauss->weight*Jdet*viscosity*thickness*(
+                                Ke->values[(4*i+3)*2*2*numnodes+4*j+1] += gauss->weight*Jdet*viscosity[2]*thickness*(
                      2.*dbasis[0*numnodes+j]*dbasis[1*numnodes+i] + dbasis[1*numnodes+j]*dbasis[0*numnodes+i]
                      )*2*pow(n+1,2)/((2*n+3)*(n+2));//K42
-                                Ke->values[(4*i+3)*2*2*numnodes+4*j+2] += gauss->weight*Jdet*viscosity*thickness*(
+                                Ke->values[(4*i+3)*2*2*numnodes+4*j+2] += gauss->weight*Jdet*viscosity[1]*thickness*(
                                                         4.*dbasis[1*numnodes+j]*dbasis[1*numnodes+i] + dbasis[0*numnodes+j]*dbasis[0*numnodes+i]
                                                         )*(n+1)/(n+2);//K43
-                                Ke->values[(4*i+3)*2*2*numnodes+4*j+3] += gauss->weight*Jdet*viscosity*thickness*(
+                                Ke->values[(4*i+3)*2*2*numnodes+4*j+3] += gauss->weight*Jdet*viscosity[2]*thickness*(
                                                         4.*dbasis[1*numnodes+j]*dbasis[1*numnodes+i] + dbasis[0*numnodes+j]*dbasis[0*numnodes+i]
                                                         )*2*pow(n+1,2)/((2*n+3)*(n+2))
                                                         + 
-                                                        gauss->weight*Jdet*viscosity*basis[j]*basis[i]*pow(n+1,2)/(thickness*(2*n+1))
+                                                        gauss->weight*Jdet*viscosity[3]*basis[j]*basis[i]*pow(n+1,2)/(thickness*(2*n+1))
                                                         ;//K44
                         }

issm/trunk-jpl/src/c/classes/Elements/Element.cpp ¶

@@ -4212 +4212 @@
 
 }/*}}}*/
-void       Element::StrainRateMLHO(IssmDouble* epsilon,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input,Input* vxshear_input,Input* vyshear_input,Input* thickness_input){/*{{{*/
-        /*Compute the 2d Blatter/MLHO Strain Rate (5 components):
+void       Element::StrainRateMLHO(IssmDouble* epsilon,IssmDouble* xyz_list,Gauss* gauss,Input* vxbase_input,Input* vybase_input,Input* vxshear_input,Input* vyshear_input,Input* thickness_input,Input* n_input,IssmDouble zeta){/*{{{*/
+        /*Compute the 2d Blatter/MLHO Strain Rate (5 components) for a given vertical coordinate (zeta):
          *
          * epsilon=[exx eyy exy exz eyz]
@@ -4221 +4221 @@
          *
          * the contribution of vz is neglected
+         *
+         * zeta = (surface - z)/thickness 
+         * i.e., zeta=0 <=> ice surface
+         * i.e., zeta=1 <=> ice base
          */
 
         /*Intermediaries*/
-        IssmDouble dvx[2]; // [dvx/dx, dvx/dy]
-        IssmDouble dvy[2]; // [dvy/dx, dvy/dy]
-        IssmDouble vxshear,vyshear,thickness;
+        IssmDouble dvxb[2];  // [dvx/dx, dvx/dy] at the base
+        IssmDouble dvyb[2];  // [dvy/dx, dvy/dy] at the base
+        IssmDouble dvxsh[2]; // [dvx/dx, dvx/dy] at the surface (only the shear)
+        IssmDouble dvysh[2]; // [dvy/dx, dvy/dy] at the surface (only the shear)
+        IssmDouble vxsh,vysh; // shear vx and shear at the surface
+        IssmDouble thickness,n;
 
         /*Check that both inputs have been found*/
-        if (!vx_input || !vy_input || !vxshear_input || !vyshear_input || !thickness_input){
-                _error_("Input missing. Here are the input pointers we have for vx: " << vx_input << ", vy: " << vy_input << ", vxshear: " << vxshear_input << ", vyshear: " << vyshear_input << ", thickness: " << thickness_input << "\n");
+        if (!vxbase_input || !vybase_input || !vxshear_input || !vyshear_input || !thickness_input || !n_input){
+                _error_("Input missing. Here are the input pointers we have for vxb: " << vxbase_input << ", vyb: " << vybase_input << ", vxshear: " << vxshear_input << ", vyshear: " << vyshear_input << ", thickness: " << thickness_input << ", n: " << n_input<< "\n");
         }
 
         /*Get strain rate assuming that epsilon has been allocated*/
-        vx_input->GetInputDerivativeValue(&dvx[0],xyz_list,gauss);
-        vy_input->GetInputDerivativeValue(&dvy[0],xyz_list,gauss);
+        vxbase_input->GetInputDerivativeValue(&dvxb[0],xyz_list,gauss);
+        vybase_input->GetInputDerivativeValue(&dvyb[0],xyz_list,gauss);
+        vxshear_input->GetInputDerivativeValue(&dvxsh[0],xyz_list,gauss);
+        vyshear_input->GetInputDerivativeValue(&dvysh[0],xyz_list,gauss);
         thickness_input->GetInputValue(&thickness,gauss); _assert_(thickness>0);
-   vxshear_input->GetInputValue(&vxshear,gauss);
-   vyshear_input->GetInputValue(&vyshear,gauss);
-        
-        epsilon[0] = dvx[0];
-        epsilon[1] = dvy[1];
-        epsilon[2] = 0.5*(dvx[1] + dvy[0]);
-        epsilon[3] = 0.5*vxshear/thickness;
-        epsilon[4] = 0.5*vyshear/thickness;
+   vxshear_input->GetInputValue(&vxsh,gauss);
+   vyshear_input->GetInputValue(&vysh,gauss);
+   n_input->GetInputValue(&n,gauss);
+
+        epsilon[0] = dvxb[0] + dvxsh[0]*( 1 - pow(zeta,n+1) );
+        epsilon[1] = dvyb[1] + dvysh[1]*( 1 - pow(zeta,n+1) );
+        epsilon[2] = 0.5*( dvxb[1] + dvxsh[1]*( 1 - pow(zeta,n+1) ) + dvyb[0] + dvysh[0]*( 1 - pow(zeta,n+1) ) ); 
+        epsilon[3] = 0.5*(vxsh/thickness)*(n+1)*pow(zeta,n);
+        epsilon[4] = 0.5*(vysh/thickness)*(n+1)*pow(zeta,n);
 
 }/*}}}*/

issm/trunk-jpl/src/c/classes/Elements/Element.h ¶

@@ -181 +181 @@
                 void               StrainRateHO(IssmDouble* epsilon,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input);
                 void               StrainRateHO2dvertical(IssmDouble* epsilon,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input);
-                void               StrainRateMLHO(IssmDouble* epsilon,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input,Input* vxshear_input,Input* vyshear_input,Input* thickness_input);
+                void               StrainRateMLHO(IssmDouble* epsilon,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input,Input* vxshear_input,Input* vyshear_input,Input* thickness_input,Input* n_input,IssmDouble zeta);
                 void               StrainRateSSA(IssmDouble* epsilon,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input);
                 void               StrainRateSSA1d(IssmDouble* epsilon,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input);

issm/trunk-jpl/src/c/classes/Materials/Material.h ¶

@@ -51 +51 @@
                 virtual void       ViscosityFS(IssmDouble* pviscosity,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input,Input* vz_input)=0;
                 virtual void       ViscosityHO(IssmDouble* pviscosity,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input)=0;
-                virtual void       ViscosityMLHO(IssmDouble* pviscosity,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input,Input* vxshear_input,Input* vyshear_input,Input* thickness_input)=0;
+                virtual void       ViscosityMLHO(IssmDouble* pviscosity,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input,Input* vxshear_input,Input* vyshear_input,Input* thickness_input,Input* n_input)=0;
                 virtual void       ViscosityL1L2(IssmDouble* pviscosity,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input,Input* surf)=0;
                 virtual void       ViscositySSA(IssmDouble* pviscosity,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input)=0;

issm/trunk-jpl/src/c/classes/Materials/Matestar.cpp ¶

@@ -567 +567 @@
         *pviscosity=GetViscosityGeneral(vx,vy,vz,&dvx[0],&dvy[0],&dvz[0],eps_eff,isdepthaveraged,gauss);
 }/*}}}*/
-void  Matestar::ViscosityMLHO(IssmDouble* pviscosity,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input,Input* vxshear_input,Input* vyshear_input,Input* thickness_input){/*{{{*/
+void  Matestar::ViscosityMLHO(IssmDouble* pviscosity,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input,Input* vxshear_input,Input* vyshear_input,Input* thickness_input,Input* n_input){/*{{{*/
         _error_("not implemented yet");
 }/*}}}*/

issm/trunk-jpl/src/c/classes/Materials/Matestar.h ¶

@@ -76 +76 @@
                 void       ViscosityFS(IssmDouble* pviscosity,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input,Input* vz_input);
                 void       ViscosityHO(IssmDouble* pviscosity,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input);
-                void       ViscosityMLHO(IssmDouble* pviscosity,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input,Input* vxshear_input,Input* vyshear_input,Input* thickness_input);
+                void       ViscosityMLHO(IssmDouble* pviscosity,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input,Input* vxshear_input,Input* vyshear_input,Input* thickness_input,Input* n_input);
                 void       ViscosityL1L2(IssmDouble* pviscosity,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input,Input* surf);
                 void       ViscositySSA(IssmDouble* pviscosity,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input);

issm/trunk-jpl/src/c/classes/Materials/Matice.cpp ¶

@@ -736 +736 @@
         *pviscosity=viscosity;
 }/*}}}*/
-void  Matice::ViscosityMLHO(IssmDouble* pviscosity,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input,Input* vxshear_input,Input* vyshear_input,Input* thickness_input){/*{{{*/
+void  Matice::ViscosityMLHO(IssmDouble* pviscosity,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vxbase_input,Input* vybase_input,Input* vxshear_input,Input* vyshear_input,Input* thickness_input,Input* n_input){/*{{{*/
 
         /*Intermediaries*/
-        IssmDouble viscosity;
-        IssmDouble thickness;
         IssmDouble epsilon[5];  /* epsilon=[exx,eyy,exy,exz,eyz]; */
-        IssmDouble eps_eff;
-
-        /* eps_eff^2 = exx^2 + eyy^2 + exy^2 + exz^2 + eyz^2 + exx*eyy */
-        element->StrainRateMLHO(&epsilon[0],xyz_list,gauss,vx_input,vy_input,vxshear_input,vyshear_input,thickness_input);
-        eps_eff = sqrt(epsilon[0]*epsilon[0] + epsilon[1]*epsilon[1] + epsilon[2]*epsilon[2] +  epsilon[3]*epsilon[3] + epsilon[4]*epsilon[4] + epsilon[0]*epsilon[1]);
-
-        /*Get viscosity*/
-        this->GetViscosityBar(&viscosity,eps_eff,gauss);
+        IssmDouble epsilon_eff;
+        IssmDouble zeta,H,n;
+        IssmDouble f[4],F[4];
+        IssmDouble mubar[4];
+        IssmDouble mu;
+        int order=5; 
+
+        for(int i=0;i<4;++i) mubar[i]=0;
+
+        GaussSeg* gauss_seg=new GaussSeg(order);
+        //IssmDouble eps_eff_averaged=0;
+        while(gauss_seg->next()){
+                
+                /*Compute zeta for gauss_seg point (0=surface, 1=base)*/
+                zeta=0.5*(gauss_seg->coord1+1); 
+
+                /* eps_eff^2 = exx^2 + eyy^2 + exy^2 + exz^2 + eyz^2 + exx*eyy (for a given zeta)*/
+                element->StrainRateMLHO(&epsilon[0],xyz_list,gauss,
+                                                vxbase_input,vybase_input,vxshear_input,vyshear_input,thickness_input,n_input,zeta);
+                epsilon_eff=sqrt(epsilon[0]*epsilon[0] + epsilon[1]*epsilon[1] + epsilon[2]*epsilon[2] 
+                                                  +  epsilon[3]*epsilon[3] + epsilon[4]*epsilon[4] + epsilon[0]*epsilon[1]);
+
+                /*Get viscosity at zeta - FIXME for now Bbar*/
+                this->GetViscosityBar(&mu,epsilon_eff,gauss);
+
+                thickness_input->GetInputValue(&H, gauss);
+      n_input->GetInputValue(&n,gauss);
+
+                /*Compute fi and Fi at zeta*/
+                f[0]=1;
+                f[1]=(1-pow(zeta,n+1));
+                f[2]=(1-pow(zeta,n+1))*(1-pow(zeta,n+1));
+                f[3]=((n+1)/H)*pow(zeta,n) * ((n+1)/H)*pow(zeta,n);
+        
+                F[0]=H;
+                F[1]=H*(n+1)/(n+2);
+                F[2]=2*H*(n+1)*(n+1)/( (2*n+3)*(n+2) );
+                F[3]=(n+1)*(n+1)/( H*(2*n+1) );
+
+                /*Sum the viscosity*/
+                mubar[0]+=(H/(2*F[0]))*gauss_seg->weight*mu*f[0];
+                mubar[1]+=(H/(2*F[1]))*gauss_seg->weight*mu*f[1];
+                mubar[2]+=(H/(2*F[2]))*gauss_seg->weight*mu*f[2];
+                mubar[3]+=(H/(2*F[3]))*gauss_seg->weight*mu*f[3];
+
+        }//while
 
         /*Assign output pointer*/
-   *pviscosity=viscosity;
+        pviscosity[0]=mubar[0];
+        pviscosity[1]=mubar[1];
+        pviscosity[2]=mubar[2];
+        pviscosity[3]=mubar[3];
+
+        /*Clean up*/
+        delete gauss_seg;
 }/*}}}*/
 void  Matice::ViscosityL1L2(IssmDouble* pviscosity,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input,Input* surface_input){/*{{{*/

issm/trunk-jpl/src/c/classes/Materials/Matice.h ¶

@@ -78 +78 @@
                 void       ViscosityFS(IssmDouble* pviscosity,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input,Input* vz_input);
                 void       ViscosityHO(IssmDouble* pviscosity,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input);
-                void       ViscosityMLHO(IssmDouble* pviscosity,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input,Input* vxshear_input,Input* vyshear_input,Input* thickness_input);
+                void       ViscosityMLHO(IssmDouble* pviscosity,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input,Input* vxshear_input,Input* vyshear_input,Input* thickness_input,Input* n_input);
                 void       ViscosityL1L2(IssmDouble* pviscosity,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input,Input* surf);
                 void       ViscositySSA(IssmDouble* pviscosity,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input);

issm/trunk-jpl/src/c/classes/Materials/Matlitho.h ¶

@@ -72 +72 @@
                 void       ViscosityFS(IssmDouble* pviscosity,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input,Input* vz_input){_error_("not supported");};
                 void       ViscosityHO(IssmDouble* pviscosity,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input){_error_("not supported");};
-                void       ViscosityMLHO(IssmDouble* pviscosity,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input,Input* vxshear_input,Input* vyshear_input,Input* thickness_input){_error_("not supported");};
+                void       ViscosityMLHO(IssmDouble* pviscosity,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input,Input* vxshear_input,Input* vyshear_input,Input* thickness_input,Input* n_input){_error_("not supported");};
                 void       ViscosityL1L2(IssmDouble* pviscosity,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input,Input* surf){_error_("not supported");};
                 void       ViscositySSA(IssmDouble* pviscosity,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input){_error_("not supported");};