Changeset 17515

Timestamp:

03/21/14 09:32:05 (11 years ago)

Author:

Mathieu Morlighem

Message:

NEW: working on FS Tetra:

Location:

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

Files:

• Tetra.cpp (modified) (3 diffs)
• Tetra.h (modified) (3 diffs)
• TetraRef.cpp (modified) (1 diff)

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

@@ -339 +339 @@
                 return false;
         }
+}
+/*}}}*/
+/*FUNCTION Tetra::HasNodeOnBed           THIS ONE{{{*/
+bool Tetra::HasNodeOnBed(){
+
+        IssmDouble values[NUMVERTICES];
+        IssmDouble sum;
+
+        /*Retrieve all inputs and parameters*/
+        GetInputListOnVertices(&values[0],MeshVertexonbedEnum);
+        sum = values[0]+values[1]+values[2]+values[3];
+
+        return sum>0.;
 }
 /*}}}*/
@@ -682 +695 @@
 }
 /*}}}*/
+/*FUNCTION Tetra::ResetFSBasalBoundaryCondition {{{*/
+void  Tetra::ResetFSBasalBoundaryCondition(void){
+
+        int numnodes = this->GetNumberOfNodes();
+
+        int          approximation;
+        IssmDouble*  vertexonbed= NULL;
+        IssmDouble   slopex,slopey,groundedice;
+        IssmDouble   xz_plane[6];
+
+        /*For FS only: we want the CS to be tangential to the bedrock*/
+        inputs->GetInputValue(&approximation,ApproximationEnum);
+        if(IsFloating() || !HasNodeOnBed() ||  approximation!=FSApproximationEnum) return;
+
+        //printf("element number %i \n",this->id);
+        /*Get inputs*/
+        Input* slopex_input=inputs->GetInput(BedSlopeXEnum); _assert_(slopex_input);
+        Input* slopey_input=inputs->GetInput(BedSlopeYEnum); _assert_(slopey_input);
+        Input* groundedicelevelset_input=inputs->GetInput(MaskGroundediceLevelsetEnum); _assert_(groundedicelevelset_input);
+        vertexonbed = xNew<IssmDouble>(numnodes);
+        this->GetInputListOnNodesVelocity(&vertexonbed[0],MeshVertexonbedEnum);
+
+        /*Loop over basal nodes and update their CS*/
+        GaussTetra* gauss = new GaussTetra();
+        for(int i=0;i<this->NumberofNodesVelocity();i++){
+
+                if(vertexonbed[i]==1){
+                        gauss->GaussNode(this->VelocityInterpolation(),i);
+
+                        slopex_input->GetInputValue(&slopex,gauss);
+                        slopey_input->GetInputValue(&slopey,gauss);
+                        groundedicelevelset_input->GetInputValue(&groundedice,gauss);
+
+                        /*New X axis          New Z axis*/
+                        xz_plane[0]=1.;       xz_plane[3]=-slopex;  
+                        xz_plane[1]=0.;       xz_plane[4]=-slopey;  
+                        xz_plane[2]=slopex;   xz_plane[5]=1.;          
+
+                        if(groundedice>0){
+                                if(this->nodes[i]->GetApproximation()==FSvelocityEnum){
+                                        this->nodes[i]->DofInSSet(2); //vz 
+                                }
+                                else _error_("Flow equation approximation"<<EnumToStringx(this->nodes[i]->GetApproximation())<<" not supported yet");
+                        }
+                        else{
+                                if(this->nodes[i]->GetApproximation()==FSvelocityEnum){
+                                        this->nodes[i]->DofInFSet(2); //vz
+                                }
+                                else _error_("Flow equation approximation"<<EnumToStringx(this->nodes[i]->GetApproximation())<<" not supported yet");
+                        }
+
+                        XZvectorsToCoordinateSystem(&this->nodes[i]->coord_system[0][0],&xz_plane[0]);
+                }
+        }
+
+        /*cleanup*/
+        xDelete<IssmDouble>(vertexonbed);
+        delete gauss;
+}
+/*}}}*/
 /*FUNCTION Tetra::SetCurrentConfiguration {{{*/
 void  Tetra::SetCurrentConfiguration(Elements* elementsin, Loads* loadsin, Nodes* nodesin, Materials* materialsin, Parameters* parametersin){
@@ -913 +986 @@
 }
 /*}}}*/
+/*FUNCTION Tetra::VelocityInterpolation{{{*/
+int Tetra::VelocityInterpolation(void){
+        return TetraRef::VelocityInterpolation();
+}
+/*}}}*/
+/*FUNCTION Tetra::PressureInterpolation{{{*/
+int Tetra::PressureInterpolation(void){
+        return TetraRef::PressureInterpolation();
+}
+/*}}}*/
 /*FUNCTION Tetra::ZeroLevelsetCoordinates{{{*/
 void Tetra::ZeroLevelsetCoordinates(IssmDouble** pxyz_zero,IssmDouble* xyz_list,int levelsetenum){

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

@@ -88 +88 @@
                 bool        HasFaceOnBed();
                 bool        HasFaceOnSurface();
+                bool        HasNodeOnBed();
                 int         Sid();
                 bool        IsOnBed();
@@ -117 +118 @@
                 IssmDouble  StabilizationParameter(IssmDouble u, IssmDouble v, IssmDouble w, IssmDouble diameter, IssmDouble kappa){_error_("not implemented yet");};
                 IssmDouble  PureIceEnthalpy(IssmDouble pressure){_error_("not implemented yet");};
-                int         PressureInterpolation(void){_error_("not implemented yet");};
+                int         PressureInterpolation(void);
                 void        ValueP1OnGauss(IssmDouble* pvalue,IssmDouble* values,Gauss* gauss){_error_("not implemented yet");};
                 void        ValueP1DerivativesOnGauss(IssmDouble* dvalue,IssmDouble* values,IssmDouble* xyz_list,Gauss* gauss){_error_("not implemented yet");};
-                int         VelocityInterpolation(void){_error_("not implemented yet");};
+                int         VelocityInterpolation(void);
                 void        GetGroundedPart(int* point1,IssmDouble* fraction1, IssmDouble* fraction2,bool* mainlyfloating){_error_("not implemented yet");};
                 IssmDouble  GetGroundedPortion(IssmDouble* xyz_list){_error_("not implemented yet");};
@@ -154 +155 @@
                 int         NodalValue(IssmDouble* pvalue, int index, int natureofdataenum){_error_("not implemented yet");};
                 void        PositiveDegreeDay(IssmDouble* pdds,IssmDouble* pds,IssmDouble signorm){_error_("not implemented yet");};
-                void        ResetFSBasalBoundaryCondition(void){_error_("not implemented yet");};
+                void        ResetFSBasalBoundaryCondition(void);
                 void        ReduceMatrices(ElementMatrix* Ke,ElementVector* pe);
                 void        SetTemporaryElementType(int element_type_in){_error_("not implemented yet");};

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

@@ -545 +545 @@
                 case P1bubblecondensedEnum: return NUMNODESP1b;
                 case P2Enum:                return NUMNODESP2;
+                case P1P1Enum:              return NUMNODESP1*2;
+                case P1P1GLSEnum:           return NUMNODESP1*2;
+                case MINIcondensedEnum:     return NUMNODESP1b+NUMNODESP1;
+                case MINIEnum:              return NUMNODESP1b+NUMNODESP1;
+                case TaylorHoodEnum:        return NUMNODESP2+NUMNODESP1;
                 default: _error_("Element type "<<EnumToStringx(this->element_type)<<" not supported yet");
         }