Changeset 16897

Timestamp:

11/22/13 15:24:52 (11 years ago)

Author:

Mathieu Morlighem

Message:

CHG: fixing thermal but still not there yet

File:

• issm/trunk-jpl/src/c/analyses/ThermalAnalysis.cpp (modified) (4 diffs)

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

@@ -231 +231 @@
                                 for(int j=0;j<numnodes;j++){
                                         Ke->values[i*numnodes+j]+=tau_parameter*D_scalar*
-                                          ((u-um)*dbasis[0*3+i]+(v-vm)*dbasis[1*3+i]+(w-wm)*dbasis[2*3+i])*((u-um)*dbasis[0*3+j]+(v-vm)*dbasis[1*3+j]+(w-wm)*dbasis[2*3+j]);
+                                          ((u-um)*dbasis[0*numnodes+i]+(v-vm)*dbasis[1*numnodes+i]+(w-wm)*dbasis[2*numnodes+i])*((u-um)*dbasis[0*numnodes+j]+(v-vm)*dbasis[1*numnodes+j]+(w-wm)*dbasis[2*numnodes+j]);
                                 }
                         }
@@ -237 +237 @@
                                 for(int i=0;i<numnodes;i++){
                                         for(int j=0;j<numnodes;j++){
-                                                Ke->values[i*numnodes+j]+=tau_parameter*D_scalar*basis[j]*((u-um)*dbasis[0*3+i]+(v-vm)*dbasis[1*3+i]+(w-wm)*dbasis[2*3+i]);
+                                                Ke->values[i*numnodes+j]+=tau_parameter*D_scalar*basis[j]*((u-um)*dbasis[0*numnodes+i]+(v-vm)*dbasis[1*numnodes+i]+(w-wm)*dbasis[2*numnodes+i]);
                                         }
                                 }
@@ -377 +377 @@
                         tau_parameter=element->StabilizationParameter(u,v,w,diameter,kappa);
 
-                        for(int i=0;i<numnodes;i++) pe->values[i]+=tau_parameter*scalar_def*(u*dbasis[0*3+i]+v*dbasis[1*3+i]+w*dbasis[2*3+i]);
+                        for(int i=0;i<numnodes;i++) pe->values[i]+=tau_parameter*scalar_def*(u*dbasis[0*numnodes+i]+v*dbasis[1*numnodes+i]+w*dbasis[2*numnodes+i]);
                         if(reCast<bool,IssmDouble>(dt)){
-                                for(int i=0;i<numnodes;i++) pe->values[i]+=tau_parameter*scalar_transient*(u*dbasis[0*3+i]+v*dbasis[1*3+i]+w*dbasis[2*3+i]);
+                                for(int i=0;i<numnodes;i++) pe->values[i]+=tau_parameter*scalar_transient*(u*dbasis[0*numnodes+i]+v*dbasis[1*numnodes+i]+w*dbasis[2*numnodes+i]);
                         }
                 }
@@ -393 +393 @@
 }/*}}}*/
 ElementVector* ThermalAnalysis::CreatePVectorSheet(Element* element){/*{{{*/
-        return NULL;
+
+        /* Geothermal flux on ice sheet base and basal friction */
+        if(!element->IsOnBed() || element->IsFloating()) return NULL;
+
+        IssmDouble  dt,Jdet,geothermalflux,vx,vy,vz;
+        IssmDouble  alpha2,scalar,basalfriction,heatflux;
+        IssmDouble *xyz_list_base = NULL;
+
+        /*Fetch number of nodes for this finite element*/
+        int numnodes = element->GetNumberOfNodes();
+
+        /*Initialize vectors*/
+        ElementVector* pe    = element->NewElementVector();
+        IssmDouble*    basis = xNew<IssmDouble>(numnodes);
+
+        /*Retrieve all inputs and parameters*/
+        element->GetVerticesCoordinatesBase(&xyz_list_base);
+        element->FindParam(&dt,TimesteppingTimeStepEnum);
+        Input* vx_input             = element->GetInput(VxEnum);                          _assert_(vx_input);
+        Input* vy_input             = element->GetInput(VyEnum);                          _assert_(vy_input);
+        Input* vz_input             = element->GetInput(VzEnum);                          _assert_(vz_input);
+        Input* geothermalflux_input = element->GetInput(BasalforcingsGeothermalfluxEnum); _assert_(geothermalflux_input);
+        IssmDouble  rho_ice             = element->GetMaterialParameter(MaterialsRhoIceEnum);
+        IssmDouble  heatcapacity        = element->GetMaterialParameter(MaterialsHeatcapacityEnum);
+
+        /*Build friction element, needed later: */
+        Friction* friction=new Friction(element,3);
+
+        /* Start  looping on the number of gaussian points: */
+        Gauss* gauss   = element->NewGaussBase(2);
+        for(int ig=gauss->begin();ig<gauss->end();ig++){
+                gauss->GaussPoint(ig);
+
+                element->JacobianDeterminantBase(&Jdet,xyz_list_base,gauss);
+                element->NodalFunctions(basis,gauss);
+
+                geothermalflux_input->GetInputValue(&geothermalflux,gauss);
+                friction->GetAlpha2(&alpha2,gauss,vx_input,vy_input,vz_input);
+                vx_input->GetInputValue(&vx,gauss);
+                vy_input->GetInputValue(&vy,gauss);
+                vz_input->GetInputValue(&vz,gauss);
+                vz = 0.;//FIXME
+                basalfriction = alpha2*(vx*vx + vy*vy + vz*vz);
+                heatflux      = (basalfriction+geothermalflux)/(rho_ice*heatcapacity);
+
+                scalar = gauss->weight*Jdet*heatflux;
+                if(dt!=0.) scalar=dt*scalar;
+
+                for(int i=0;i<numnodes;i++) pe->values[i]+=scalar*basis[i];
+        }
+
+        /*Clean up and return*/
+        delete gauss;
+        delete friction;
+        xDelete<IssmDouble>(basis);
+        xDelete<IssmDouble>(xyz_list_base);
+        return pe;
 }/*}}}*/
 ElementVector* ThermalAnalysis::CreatePVectorShelf(Element* element){/*{{{*/