Changeset 16866

Timestamp:

11/21/13 15:45:51 (11 years ago)

Author:

Mathieu Morlighem

Message:

CHG: added FS vertical

Location:

issm/trunk-jpl/src/c

Files:

• analyses/StressbalanceAnalysis.cpp (modified) (1 diff)
• classes/Elements/Element.h (modified) (1 diff)
• classes/Elements/Penta.h (modified) (1 diff)
• classes/Elements/Seg.h (modified) (1 diff)
• classes/Elements/Tria.cpp (modified) (5 diffs)
• classes/Elements/Tria.h (modified) (4 diffs)

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

@@ -2228 +2228 @@
                 element->NodalFunctionsVelocity(vbasis,gauss);
                 surface_input->GetInputValue(&surface,gauss);
-                z=element->GetZcoord(gauss);
+                if(dim==3) z=element->GetZcoord(gauss);
+                else       z=element->GetYcoord(gauss);
                 pressure = rho_water*gravity*min(0.,z);//0 if the gaussian point is above water level
 

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

@@ -120 +120 @@
                 virtual void   GetVerticesCoordinatesTop(IssmDouble** xyz_list)=0;
                 virtual IssmDouble GetZcoord(Gauss* gauss)=0;
+                virtual IssmDouble GetYcoord(Gauss* gauss)=0;
                 virtual void   GetMaterialInputValue(IssmDouble* pvalue,Node* node,int enumtype)=0;
 

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

@@ -102 +102 @@
                 void   GetSolutionFromInputsOneDof(Vector<IssmDouble>* solution,int enum_type);
                 IssmDouble GetZcoord(Gauss* gauss);
+                IssmDouble GetYcoord(Gauss* gauss){_error_("Not implemented");};
                 void   GetVectorFromInputs(Vector<IssmDouble>* vector,int name_enum);
                 void   GetVerticesCoordinates(IssmDouble** pxyz_list);

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

@@ -151 +151 @@
                 void        GetInputValue(IssmDouble* pvalue,int enum_type){_error_("not implemented yet");};
                 void        GetMaterialInputValue(IssmDouble* pvalue,Node* node,int enumtype){_error_("not implemented yet");};
+                IssmDouble  GetYcoord(Gauss* gauss){_error_("Not implemented");};
                 IssmDouble  GetZcoord(Gauss* gauss){_error_("not implemented yet");};
                 Gauss*      NewGauss(void){_error_("not implemented yet");};

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

@@ -1190 +1190 @@
 
 }/*}}}*/
+/*FUNCTION Tria::GetVerticesCoordinatesBase(IssmDouble** pxyz_list){{{*/
+void Tria::GetVerticesCoordinatesBase(IssmDouble** pxyz_list){
+
+        int        indices[2];
+        IssmDouble xyz_list[NUMVERTICES][3];
+
+        /*Element XYZ list*/
+        ::GetVerticesCoordinates(&xyz_list[0][0],this->vertices,NUMVERTICES);
+
+        /*Allocate Output*/
+        IssmDouble* xyz_list_edge = xNew<IssmDouble>(2*3);
+        this->EdgeOnBedIndices(&indices[0],&indices[1]);
+        for(int i=0;i<2;i++) for(int j=0;j<2;j++) xyz_list_edge[i*3+j]=xyz_list[indices[i]][j];
+
+        /*Assign output pointer*/
+        *pxyz_list = xyz_list_edge;
+
+}/*}}}*/
 /*FUNCTION Tria::NormalSection{{{*/
 void Tria::NormalSection(IssmDouble* normal,IssmDouble* xyz_list){
@@ -1577 +1595 @@
         /*We found the enum.  Use its values to fill into the vector, using the vertices ids: */
         input->GetVectorFromInputs(vector,&vertexpidlist[0]);
+}
+/*}}}*/
+/*FUNCTION Tria::GetYcoord {{{*/
+IssmDouble Tria::GetYcoord(Gauss* 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];
+        int backup = this->element_type;
+        this->element_type = P1Enum;
+        TriaRef::GetInputValue(&y,y_list,gauss);
+        this->element_type = backup;
+
+        return y;
 }
 /*}}}*/
@@ -2044 +2079 @@
 }
 /*}}}*/
+/*FUNCTION Tria::JacobianDeterminantBase{{{*/
+void Tria::JacobianDeterminantBase(IssmDouble* pJdet,IssmDouble* xyz_list_base,Gauss* gauss){
+
+        _assert_(gauss->Enum()==GaussTriaEnum);
+        this->GetSegmentJacobianDeterminant(pJdet,xyz_list_base,(GaussTria*)gauss);
+
+}
+/*}}}*/
 /*FUNCTION Tria::JacobianDeterminantSurface{{{*/
 void Tria::JacobianDeterminantSurface(IssmDouble* pJdet,IssmDouble* xyz_list,Gauss* gauss){
@@ -2263 +2306 @@
         GetAreaCoordinates(&area_coordinates[0][0],xyz_list_front,xyz_list,4);
         return new GaussTria(area_coordinates,order_vert);
+}
+/*}}}*/
+/*FUNCTION Tria::NewGaussBase(int order){{{*/
+Gauss* Tria::NewGaussBase(int order){
+
+        int indices[2];
+        this->EdgeOnBedIndices(&indices[0],&indices[1]);
+        return new GaussTria(indices[0],indices[1],order);
 }
 /*}}}*/
@@ -2360 +2411 @@
         if(found)*pvalue=value;
         return found;
+}
+/*}}}*/
+/*FUNCTION Tria::NormalBase {{{*/
+void Tria::NormalBase(IssmDouble* bed_normal,IssmDouble* xyz_list){
+
+        /*Build unit outward pointing vector*/
+        IssmDouble vector[2];
+        IssmDouble norm;
+
+        vector[0]=xyz_list[1*3+0] - xyz_list[0*3+0];
+        vector[1]=xyz_list[1*3+1] - xyz_list[0*3+1];
+
+        norm=sqrt(vector[0]*vector[0] + vector[1]*vector[1]);
+
+        bed_normal[0]= + vector[1]/norm;
+        bed_normal[1]= - vector[0]/norm;
 }
 /*}}}*/

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

@@ -117 +117 @@
                 void        GetVectorFromInputs(Vector<IssmDouble>* vector, int name_enum);
                 void        GetVerticesCoordinates(IssmDouble** pxyz_list);
-                void        GetVerticesCoordinatesBase(IssmDouble** pxyz_list){_error_("not implemented yet");};
+                void        GetVerticesCoordinatesBase(IssmDouble** pxyz_list);
                 void        GetVerticesCoordinatesTop(IssmDouble** pxyz_list){_error_("not implemented yet");};
                 void        InputCreate(IssmDouble* vector,IoModel* iomodel,int M,int N,int vector_type,int vector_enum,int code);
@@ -266 +266 @@
                 void           GetGroundedPart(int* point1,IssmDouble* fraction1, IssmDouble* fraction2,bool* mainlyfloating);
                 IssmDouble     GetGroundedPortion(IssmDouble* xyz_list);
+                IssmDouble     GetYcoord(Gauss* gauss);
                 IssmDouble     GetZcoord(Gauss* gauss){_error_("not implemented");};
                 void           NormalSection(IssmDouble* normal,IssmDouble* xyz_list);
                 void           NormalTop(IssmDouble* normal,IssmDouble* xyz_list){_error_("not implemented yet");};
-                void           NormalBase(IssmDouble* normal,IssmDouble* xyz_list){_error_("not implemented yet");};
+                void           NormalBase(IssmDouble* normal,IssmDouble* xyz_list);
                 IssmDouble     GetMaterialParameter(int enum_in);
                 Input*         GetInput(int inputenum);
@@ -291 +292 @@
                 void           JacobianDeterminantLine(IssmDouble* Jdet, IssmDouble* xyz_list,Gauss* gauss){_error_("not implemented yet");};
                 void           JacobianDeterminantSurface(IssmDouble*  pJdet, IssmDouble* xyz_list,Gauss* gauss);
-                void           JacobianDeterminantBase(IssmDouble* pJdet,IssmDouble* xyz_list_base,Gauss* gauss){_error_("not implemented yet");};
+                void           JacobianDeterminantBase(IssmDouble* pJdet,IssmDouble* xyz_list_base,Gauss* gauss);
                 void           JacobianDeterminantTop(IssmDouble* pJdet,IssmDouble* xyz_list_base,Gauss* gauss){_error_("not implemented yet");};
                 IssmDouble     MinEdgeLength(IssmDouble* xyz_list){_error_("not implemented yet");};
@@ -298 +299 @@
       Gauss*         NewGauss(IssmDouble* xyz_list, IssmDouble* xyz_list_front,int order);
       Gauss*         NewGauss(IssmDouble* xyz_list, IssmDouble* xyz_list_front,int order_horiz,int order_vert);
-                Gauss*         NewGaussBase(int order){_error_("not implemented yet");};
+                Gauss*         NewGaussBase(int order);
                 Gauss*         NewGaussLine(int vertex1,int vertex2,int order){_error_("not implemented yet");};
                 Gauss*         NewGaussTop(int order){_error_("not implemented yet");};