Changeset 17949

Timestamp:

05/06/14 13:45:41 (11 years ago)

Author:

hongjuy

Message:

ADD: shelf dampening force for basal melting

File:

• issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp (modified) (2 diffs)

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

@@ -212 +212 @@
         iomodel->FetchDataToInput(elements,MaterialsRheologyBEnum);
         iomodel->FetchDataToInput(elements,MaterialsRheologyNEnum);
+        iomodel->FetchDataToInput(elements,BasalforcingsMeltingRateEnum);
         iomodel->FetchDataToInput(elements,VxEnum,0.);
         if(dakota_analysis)elements->InputDuplicate(VxEnum,QmuVxEnum);
@@ -3597 +3598 @@
         element->TransformLoadVectorCoord(pe,cs_list);
 
+        /* shelf dampening*/
+        int shelf_dampening;
+        element->FindParam(&shelf_dampening,StressbalanceShelfDampeningEnum);
+        if(shelf_dampening) {
+                Input*      mb_input=element->GetInput(BasalforcingsMeltingRateEnum); _assert_(mb_input);
+                IssmDouble dt,mb,normal_b;
+                element->FindParam(&dt,TimesteppingTimeStepEnum);
+                for(int ig=gauss->begin();ig<gauss->end();ig++){
+                        gauss->GaussPoint(ig);
+                        element->JacobianDeterminantBase(&Jdet,xyz_list_base,gauss);
+                        element->NodalFunctionsVelocity(vbasis,gauss);
+                        element->NormalBase(&normal[0],xyz_list_base);
+                        if (dim==2) normal_b=normal[1];
+                        else if (dim==3) normal_b=sqrt(normal[0]*normal[0]+normal[1]*normal[1]);
+                        mb_input->GetInputValue(&mb, gauss);
+                        for(i=0;i<vnumnodes;i++){
+                                pe->values[i*dim+1] += dt*rho_water*gravity*mb*gauss->weight*Jdet*vbasis[i]*normal_b;
+                        }
+                }
+        }
+
         /*Clean up and return*/
         delete gauss;