Changeset 16685 for issm/trunk-jpl/src/c/classes/Elements/Penta.cpp

Timestamp:

11/08/13 11:50:50 (11 years ago)

Author:

bdef

Message:

Adding stuff to compute Epl slopes only on active EPL

File:

• issm/trunk-jpl/src/c/classes/Elements/Penta.cpp (modified) (7 diffs)

issm/trunk-jpl/src/c/classes/Elements/Penta.cpp

@@ -481 +481 @@
           case HydrologyDCEfficientAnalysisEnum:
                         return CreateKMatrixHydrologyDCEfficient();
+                        break;
+          case L2ProjectionEPLAnalysisEnum:
+                        return CreateEPLDomainMassMatrix();
                         break;
                 #endif
@@ -697 +700 @@
                         return CreatePVectorHydrologyDCEfficient();
                         break;
+          case L2ProjectionEPLAnalysisEnum:
+                        return CreatePVectorL2ProjectionEPL();
+                        break;
                 #endif
                 default:
@@ -2258 +2264 @@
         case HydrologyDCEfficientAnalysisEnum:
                 InputUpdateFromSolutionOneDofCollapsed(solution,EplHeadEnum);
+                break;
+        case L2ProjectionEPLAnalysisEnum:
+                this->parameters->FindParam(&inputenum,InputToL2ProjectEnum);
+                InputUpdateFromSolutionOneDofCollapsed(solution,inputenum);
                 break;
         #endif
@@ -10622 +10632 @@
 }
 /*}}}*/
+/*FUNCTION Penta::CreateEPLDomainMassMatrix {{{*/
+ElementMatrix* Penta::CreateEPLDomainMassMatrix(void){
+
+        if (!IsOnBed()) return NULL;
+
+        Tria* tria=(Tria*)SpawnTria(0); //lower face is 0, upper face is 1.
+        ElementMatrix* Ke=tria->CreateEPLDomainMassMatrix();
+        delete tria->material; delete tria;
+
+        /*clean up and return*/
+        return Ke;
+}
+/*}}}*/
 /*FUNCTION Penta::CreatePVectorHydrologyDCInefficient {{{*/
 ElementVector* Penta::CreatePVectorHydrologyDCInefficient(void){
@@ -10644 +10667 @@
         Tria* tria=(Tria*)SpawnTria(0); //lower face is 0, upper face is 1.
         ElementVector* pe=tria->CreatePVectorHydrologyDCEfficient();
+        delete tria->material; delete tria;
+
+        /*Clean up and return*/
+        return pe;
+}
+/*}}}*/
+/*FUNCTION Penta::CreatePVectorL@ProjectionEPL {{{*/
+ElementVector* Penta::CreatePVectorL2ProjectionEPL(void){
+
+        if (!IsOnBed()) return NULL;
+
+        /*Call Tria function*/
+        Tria* tria=(Tria*)SpawnTria(0); //lower face is 0, upper face is 1.
+        ElementVector* pe=tria->CreatePVectorL2ProjectionEPL();
         delete tria->material; delete tria;
 
@@ -10775 +10812 @@
                                 /*Compute first the effective pressure in the EPL*/
                                 EPL_N=gravity*((rho_ice*ice_thickness[i])-(rho_water*(eplhead[i]-bed[i])));
-                                
+                                if(EPL_N<0.0)EPL_N=0.0;
                                 /*Get then the gradient of EPL heads*/
                                 EPLgrad = epl_slopeX[i]+epl_slopeY[i];
                                 
                                 /*And proceed to the real thing*/
-                                thickness[i] = old_thickness[i]*(1-((rho_water*gravity*dt)/(rho_ice* latentheat))*epl_conductivity*pow(EPLgrad,2.0)-2.0*(A*dt/(pow(n,n)))*EPL_N);
+                                thickness[i] = old_thickness[i]*(1+((rho_water*gravity*dt)/(rho_ice*latentheat))*epl_conductivity*pow(EPLgrad,2.0)-2.0*(A*dt/(pow(n,n)))*(pow(EPL_N,n)));
                                 thickness[i+numdof2d]=thickness[i];
+                                printf("N, %e - thick term2, %e \n",EPL_N,-2.0*(A*dt/(pow(n,n)))*(pow(EPL_N,n)));
+                                printf("old thick, %e - thick , %e \n",old_thickness[i],thickness[i]);
                         }
                 }
@@ -10866 +10905 @@
         /*Free ressources:*/
         xDelete<int>(doflist);
+}
+/*}}}*/
+/*FUNCTION Penta::UpdateConstraintsL2ProjectionEPL{{{*/
+void  Penta::UpdateConstraintsL2ProjectionEPL(void){
+
+        IssmDouble activeEpl[NUMVERTICES];
+
+
+        if(!IsOnBed()){
+                for(int i=0;i<this->NumberofNodes();i++)this->nodes[i]->Deactivate();
+        }
+        else{
+                GetInputListOnVertices(&activeEpl[0],HydrologydcMaskEplactiveEnum);
+                for(int i=0;i<3;i++){
+                        if(!activeEpl[i])this->nodes[i]->Deactivate();
+                }
+        }
 }
 /*}}}*/