Changeset 18545

Timestamp:

09/26/14 11:19:18 (10 years ago)

Author:

Mathieu Morlighem

Message:

CHG: last change before testing

File:

• issm/trunk-jpl/src/c/analyses/Balancethickness2Analysis.cpp (modified) (7 diffs)

issm/trunk-jpl/src/c/analyses/Balancethickness2Analysis.cpp

@@ -30 +30 @@
         iomodel->FetchDataToInput(elements,BaseEnum);
         iomodel->FetchDataToInput(elements,MaskIceLevelsetEnum);
+        iomodel->FetchDataToInput(elements,FrictionCoefficientEnum);
         iomodel->FetchDataToInput(elements,BasalforcingsGroundediceMeltingRateEnum);
         iomodel->FetchDataToInput(elements,SurfaceforcingsMassBalanceEnum);
@@ -105 +106 @@
 
         /*Intermediaries */
-        IssmDouble  dhdt,mb,ms,dD[2],db[2],Jdet;
+        IssmDouble  dhdt,mb,ms,D,db[2],Jdet;
         IssmDouble* xyz_list = NULL;
 
@@ -114 +115 @@
         ElementVector* pe    = element->NewElementVector();
         IssmDouble*    basis = xNew<IssmDouble>(numnodes);
+        IssmDouble*   dbasis = xNew<IssmDouble>(2*numnodes);
 
         /*Retrieve all inputs and parameters*/
@@ -130 +132 @@
                 element->JacobianDeterminant(&Jdet,xyz_list,gauss);
                 element->NodalFunctions(basis,gauss);
+                element->NodalFunctionsDerivatives(dbasis,xyz_list,gauss);
 
                 ms_input->GetInputValue(&ms,gauss);
                 mb_input->GetInputValue(&mb,gauss);
                 dhdt_input->GetInputValue(&dhdt,gauss);
+                D_input->GetInputValue(&D,gauss);
                 bed_input->GetInputDerivativeValue(&db[0],xyz_list,gauss);
-                D_input->GetInputDerivativeValue(&dD[0],xyz_list,gauss);
-                db[0]=0.;
-                db[1]=0.;
-
-                /*Since grad(b) is constant div(D grad(b) ) = grad(D).grad(b)*/
-                for(int i=0;i<numnodes;i++) pe->values[i]+=Jdet*gauss->weight*(ms-mb-dhdt +dD[0]*db[0]+dD[1]*db[1])*basis[i];
+                //db[0]=0.;
+                //db[1]=0.;
+
+                for(int i=0;i<numnodes;i++) pe->values[i]+=Jdet*gauss->weight*(
+                                        (ms-mb-dhdt)*basis[i]
+                                        -D*db[0]*dbasis[0*numnodes+i] - D*db[1]*dbasis[1*numnodes+i]
+                                        );
         }
 
@@ -167 +172 @@
 
         /*Intermediaries */
-        IssmDouble       omega,h,mu0,ds[2],Cmu,B;
-        const int        n = 3.;
+        IssmDouble       Gamma,h,mu0,ds[2],Cmu,B,k;
+        const int        n = 3;
         const IssmDouble Hstar = 500.;
         const IssmDouble Lstar = 500.e+3;
@@ -183 +188 @@
         Input* surface_input    = element->GetInput(SurfaceEnum);            _assert_(surface_input);
         Input* B_input          = element->GetInput(MaterialsRheologyBEnum); _assert_(B_input);
+        Input* k_input          = element->GetInput(FrictionCoefficientEnum);_assert_(k_input);
         IssmDouble rhog = element->GetMaterialParameter(MaterialsRhoIceEnum)*element->GetMaterialParameter(ConstantsGEnum);
 
@@ -191 +197 @@
                 
                 B_input->GetInputValue(&B,gauss);
+                k_input->GetInputValue(&k,gauss);
                 thickness_input->GetInputValue(&h,gauss);
                 surface_input->GetInputDerivativeValue(&ds[0],xyz_list,gauss);
 
-                mu0   = pow(2.,(1.-3*n)/(2.*n)) * B;
-                omega = pow(rhog,3) * pow(Hstar,2*(n+1)) * pow(Hstar/Lstar,2*(n+1)) * 1./pow(mu0,n);
-                Cmu   = 0.;
-
-                D[i] = omega*(ds[0]*ds[0]+ds[1]*ds[1])*pow(h,4)*(1./5.*h - Cmu);
+                mu0   = pow(2.,(1-3*n)/(2.*n)) * B;
+                Gamma = pow(rhog,n) * pow(Hstar,2*(n+1)) * pow(Hstar/Lstar,2*(n+1)) * 1./pow(mu0,n);
+                Cmu   = mu0/Hstar/(k*k);
+
+                D[i] = Gamma*pow(ds[0]*ds[0]+ds[1]*ds[1],(n-1)/2)*pow(h,n+1)*(1./(n+2.)*h - Cmu);
         }