Changeset 16364

Timestamp:

10/10/13 10:00:54 (11 years ago)

Author:

Mathieu Morlighem

Message:

NEW: adding front to FS

Location:

issm/trunk-jpl/src/c/classes/Elements

Files:

• Tria.cpp (modified) (3 diffs)
• Tria.h (modified) (1 diff)
• TriaRef.cpp (modified) (4 diffs)
• TriaRef.h (modified) (2 diffs)

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

@@ -3456 +3456 @@
         if(!isfront) return NULL;
 
-        _error_("STOP");
+        /*Intermediaries*/
+        IssmDouble  rho_ice,rho_water,gravity,y_g;
+        IssmDouble  Jdet,water_pressure,air_pressure,pressure;
+        IssmDouble  xyz_list_front[2][3];
+        IssmDouble  area_coordinates[2][3];
+        IssmDouble  normal[2];
+        GaussTria*  gauss = NULL;
+
+        /*Fetch number of nodes and dof for this finite element*/
+        int vnumnodes = this->NumberofNodesVelocity();
+        int pnumnodes = this->NumberofNodesPressure();
+
+        /*Prepare coordinate system list*/
+        int* cs_list = xNew<int>(vnumnodes+pnumnodes);
+        for(i=0;i<vnumnodes;i++) cs_list[i]           = XYEnum;
+        for(i=0;i<pnumnodes;i++) cs_list[vnumnodes+i] = PressureEnum;
+
+        /*Initialize Element matrix and vectors*/
+        GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
+        ElementVector* pe     = new ElementVector(nodes,vnumnodes+pnumnodes,this->parameters,FSvelocityEnum);
+        IssmDouble*    vbasis = xNew<IssmDouble>(vnumnodes);
+
+        /*Retrieve all inputs and parameters*/
+        rho_water = matpar->GetRhoWater();
+        gravity   = matpar->GetG();
+        GetZeroLevelsetCoordinates(&xyz_list_front[0][0],xyz_list,MaskIceLevelsetEnum);
+        GetAreaCoordinates(&area_coordinates[0][0],xyz_list_front,xyz_list,2);
+        GetSegmentNormal(&normal[0],xyz_list_front);
+
+        /*Start looping on Gaussian points*/
+        IssmDouble ymax=max(xyz_list_front[0][1],xyz_list_front[1][1]);
+        IssmDouble ymin=min(xyz_list_front[0][1],xyz_list_front[1][1]);
+        if(ymax>0. && ymin<0.) gauss=new GaussTria(area_coordinates,30); //refined in vertical because of the sea level discontinuity
+        else                   gauss=new GaussTria(area_coordinates,3);
+
+        for(int ig=gauss->begin();ig<gauss->end();ig++){
+
+                gauss->GaussPoint(ig);
+                y_g=GetYcoord(gauss);
+                GetNodalFunctionsVelocity(vbasis,gauss);
+                GetSegmentJacobianDeterminant(&Jdet,&xyz_list_front[0][0],gauss);
+
+                water_pressure = rho_water*gravity*min(0.,y_g); //0 if the gaussian point is above water level
+                air_pressure   = 0.;
+                pressure       = water_pressure + air_pressure;
+
+                for(i=0;i<vnumnodes;i++){
+                        pe->values[2*i+0]+= pressure*Jdet*gauss->weight*normal[0]*vbasis[i];
+                        pe->values[2*i+1]+= pressure*Jdet*gauss->weight*normal[1]*vbasis[i];
+                }
+        }
+
+
+        /*Transform coordinate system*/
+        TransformLoadVectorCoord(pe,nodes,(vnumnodes+pnumnodes),cs_list);
+
+        /*Clean up and return*/
+        xDelete<int>(cs_list);
+        xDelete<IssmDouble>(vbasis);
+        delete gauss;
+        return pe;
+
 }
 /*}}}*/
@@ -3579 +3640 @@
         }
 
+        pe->CheckConsistency();
         /*Transform coordinate system*/
         TransformLoadVectorCoord(pe,nodes,(vnumnodes+pnumnodes),cs_list);
+        pe->CheckConsistency();
 
         /*Clean up and return*/
@@ -3985 +4048 @@
         delete gauss;
         xDelete<int>(doflist);
+}
+/*}}}*/
+/*FUNCTION Tria::GetYcoord {{{*/
+IssmDouble Tria::GetYcoord(GaussTria* gauss){
+
+        IssmDouble y;
+        IssmDouble xyz_list[NUMVERTICES][3];
+        IssmDouble y_list[NUMVERTICES];
+
+        GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
+        for(int i=0;i<NUMVERTICES;i++) y_list[i]=xyz_list[i][1];
+        TriaRef::GetInputValue(&y,&y_list[0],gauss,P1Enum);
+
+        return y;
 }
 /*}}}*/

issm/trunk-jpl/src/c/classes/Elements/Tria.h

@@ -248 +248 @@
                 ElementVector* CreatePVectorStressbalanceFSShelf(void);
                 ElementMatrix* CreateJacobianStressbalanceSSA(void);
-                void      GetSolutionFromInputsStressbalanceFS(Vector<IssmDouble>* solution);
-                void      GetSolutionFromInputsStressbalanceHoriz(Vector<IssmDouble>* solution);
-                void      GetSolutionFromInputsStressbalanceSIA(Vector<IssmDouble>* solution);
-                void      InputUpdateFromSolutionStressbalanceHoriz( IssmDouble* solution);
-                void      InputUpdateFromSolutionStressbalanceFS( IssmDouble* solution);
-                void      InputUpdateFromSolutionStressbalanceSIA( IssmDouble* solution);
+                void             GetSolutionFromInputsStressbalanceFS(Vector<IssmDouble>* solution);
+                void             GetSolutionFromInputsStressbalanceHoriz(Vector<IssmDouble>* solution);
+                void             GetSolutionFromInputsStressbalanceSIA(Vector<IssmDouble>* solution);
+                IssmDouble     GetYcoord(GaussTria* gauss);
+                void             InputUpdateFromSolutionStressbalanceHoriz( IssmDouble* solution);
+                void             InputUpdateFromSolutionStressbalanceFS( IssmDouble* solution);
+                void             InputUpdateFromSolutionStressbalanceSIA( IssmDouble* solution);
                 #endif
 

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

@@ -533 +533 @@
 }
 /*}}}*/
-/*FUNCTION TriaRef::GetNodalFunctions{{{*/
+/*FUNCTION TriaRef::GetNodalFunctions(IssmDouble* basis,GaussTria* gauss){{{*/
 void TriaRef::GetNodalFunctions(IssmDouble* basis,GaussTria* gauss){
         /*This routine returns the values of the nodal functions  at the gaussian point.*/
 
         _assert_(basis);
-
-        switch(this->element_type){
+        GetNodalFunctions(basis,gauss,this->element_type);
+
+}
+/*}}}*/
+/*FUNCTION TriaRef::GetNodalFunctions(IssmDouble* basis,GaussTria* gauss,int interpolation){{{*/
+void TriaRef::GetNodalFunctions(IssmDouble* basis,GaussTria* gauss,int interpolation){
+        /*This routine returns the values of the nodal functions  at the gaussian point.*/
+
+        _assert_(basis);
+
+        switch(interpolation){
                 case P1Enum: case P1DGEnum:
                         basis[0]=gauss->coord1;
@@ -867 +876 @@
 }
 /*}}}*/
-/*FUNCTION TriaRef::GetInputValue{{{*/
+/*FUNCTION TriaRef::GetInputValue(IssmDouble* p, IssmDouble* plist, GaussTria* gauss){{{*/
 void TriaRef::GetInputValue(IssmDouble* p, IssmDouble* plist, GaussTria* gauss){
+
+        GetInputValue(p,plist,gauss,this->element_type);
+}
+/*}}}*/
+/*FUNCTION TriaRef::GetInputValue(IssmDouble* p, IssmDouble* plist, GaussTria* gauss,int interpolation){{{*/
+void TriaRef::GetInputValue(IssmDouble* p, IssmDouble* plist, GaussTria* gauss,int interpolation){
 
         /*Output*/
@@ -874 +889 @@
 
         /*Fetch number of nodes for this finite element*/
-        int numnodes = this->NumberofNodes();
+        int numnodes = this->NumberofNodes(interpolation);
 
         /*Get nodal functions*/
         IssmDouble* basis=xNew<IssmDouble>(numnodes);
-        GetNodalFunctions(basis, gauss);
+        GetNodalFunctions(basis, gauss,interpolation);
 
         /*Calculate parameter for this Gauss point*/
@@ -888 +903 @@
 }
 /*}}}*/
-/*FUNCTION TriaRef::NumberofNodes{{{*/
+/*FUNCTION TriaRef::NumberofNodes(){{{*/
 int TriaRef::NumberofNodes(void){
 
-        switch(this->element_type){
+        return this->NumberofNodes(this->element_type);
+}
+/*}}}*/
+/*FUNCTION TriaRef::NumberofNodes(int interpolation){{{*/
+int TriaRef::NumberofNodes(int interpolation){
+
+        switch(interpolation){
                 case P1Enum:                return NUMNODESP1;
                 case P1DGEnum:              return NUMNODESP1;

issm/trunk-jpl/src/c/classes/Elements/TriaRef.h

@@ -38 +38 @@
                 void GetJacobianInvert(IssmDouble*  Jinv, IssmDouble* xyz_list,GaussTria* gauss);
                 void GetNodalFunctions(IssmDouble* basis,GaussTria* gauss);
+                void GetNodalFunctions(IssmDouble* basis,GaussTria* gauss,int interpolation);
                 void GetNodalFunctionsVelocity(IssmDouble* basis, GaussTria* gauss);
                 void GetNodalFunctionsPressure(IssmDouble* basis, GaussTria* gauss);
@@ -48 +49 @@
                 void GetNodalFunctionsDerivativesReference(IssmDouble* dbasis,GaussTria* gauss);
                 void GetInputValue(IssmDouble* pp, IssmDouble* plist, GaussTria* gauss);
+                void GetInputValue(IssmDouble* pp, IssmDouble* plist, GaussTria* gauss,int interpolation);
                 void GetInputDerivativeValue(IssmDouble* pp, IssmDouble* plist,IssmDouble* xyz_list, GaussTria* gauss);
 
                 int  NumberofNodes(void);
+                int  NumberofNodes(int interpolation);
                 int  NumberofNodesVelocity(void);
                 int  NumberofNodesPressure(void);