Changeset 16387

Timestamp:

10/11/13 11:48:46 (11 years ago)

Author:

Mathieu Morlighem

Message:

CHG: Fixing local coordinate system for flowband model

Location:

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

Files:

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

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

@@ -2182 +2182 @@
 }
 /*}}}*/
+/*FUNCTION Tria::EdgeOnBedIndex{{{*/
+int Tria::EdgeOnBedIndex(void){
+
+        IssmDouble values[NUMVERTICES];
+        int        indices[3][2] = {{1,2},{2,0},{0,1}};
+
+        /*Retrieve all inputs and parameters*/
+        GetInputListOnVertices(&values[0],MeshVertexonbedEnum);
+
+        for(int i=0;i<3;i++){
+                if(values[indices[i][0]] == 1. && values[indices[i][1]] == 1.){
+                        return i;
+                }
+        }
+
+        _printf_("list of vertices on bed: "<<values[0]<<" "<<values[1]<<" "<<values[2]);
+        _error_("Could not find 2 vertices on bed");
+}
+/*}}}*/
 /*FUNCTION Tria::IsFloating {{{*/
 bool   Tria::IsFloating(){
@@ -2408 +2427 @@
 
         /*clean-up*/
+        delete gauss;
+}
+/*}}}*/
+/*FUNCTION Tria::ResetCoordinateSystem{{{*/
+void  Tria::ResetCoordinateSystem(void){
+
+        int        approximation;
+        int        numindices;
+        int       *indices = NULL;
+        IssmDouble slope;
+        IssmDouble xz_plane[6];
+
+        /*For FS only: we want the CS to be tangential to the bedrock*/
+        inputs->GetInputValue(&approximation,ApproximationEnum);
+        if(IsFloating() || !HasEdgeOnBed() || approximation!=FSApproximationEnum) return;
+
+        /*Get number of nodes for velocity only and base*/
+        int index = this->EdgeOnBedIndex();
+        NodeOnEdgeIndices(&numindices,&indices,index,this->VelocityInterpolation());
+
+        /*Get inputs*/
+        Input* slope_input=inputs->GetInput(BedSlopeXEnum); _assert_(slope_input);
+
+        /*Loop over basal nodes and update their CS*/
+        GaussTria* gauss = new GaussTria();
+        for(int i=0;i<numindices;i++){//FIXME
+
+
+                gauss->GaussNode(this->VelocityInterpolation(),indices[i]);
+
+                slope_input->GetInputValue(&slope,gauss);
+
+                IssmDouble theta = atan(slope);
+
+                /*New X axis                  New Z axis*/
+                xz_plane[0]=cos(theta);       xz_plane[3]=0.;  
+                xz_plane[1]=sin(theta);       xz_plane[4]=0.;  
+                xz_plane[2]=0.;               xz_plane[5]=1.;          
+
+                XZvectorsToCoordinateSystem(&this->nodes[indices[i]]->coord_system[0][0],&xz_plane[0]);
+        }
+
+        /*cleanup*/
+        xDelete<int>(indices);
         delete gauss;
 }

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

@@ -87 +87 @@
                 bool        HasEdgeOnBed();
                 void        EdgeOnBedIndices(int* pindex1,int* pindex);
+                int         EdgeOnBedIndex();
                 bool        IsFloating();
                 bool        IsNodeOnShelfFromFlags(IssmDouble* flags);
@@ -106 +107 @@
                 void        PatchFill(int* pcount, Patch* patch);
                 void        PatchSize(int* pnumrows, int* pnumvertices,int* pnumnodes);
-                void        ResetCoordinateSystem(void){_error_("not implemented yet");};
+                void        ResetCoordinateSystem(void);
                 void          SmbGradients();
                 IssmDouble  SurfaceArea(void);

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

@@ -1038 +1038 @@
 }
 /*}}}*/
+/*FUNCTION TriaRef::NodeOnEdgeIndices{{{*/
+void TriaRef::NodeOnEdgeIndices(int* pnumindices,int** pindices,int index,int finiteelement){
+
+        /*Output*/
+        int  numindices;
+        int* indices = NULL;
+
+        switch(finiteelement){
+                case P1Enum: case P1DGEnum: case P1bubbleEnum: case P1bubblecondensedEnum:
+                        numindices = 2;
+                        indices    = xNew<int>(numindices);
+                        switch(index){
+                                case 0:
+                                        indices[0] = 1;
+                                        indices[1] = 2;
+                                        break;
+                                case 1:
+                                        indices[0] = 2;
+                                        indices[1] = 0;
+                                        break;
+                                case 2:
+                                        indices[0] = 0;
+                                        indices[1] = 1;
+                                        break;
+                                default:
+                                        _error_("Edge index provided ("<<index<<") is not between 0 and 2");
+                        }
+                        break;
+                case P2Enum:
+                        numindices = 3;
+                        indices    = xNew<int>(numindices);
+                        switch(index){
+                                case 0:
+                                        indices[0] = 1;
+                                        indices[1] = 2;
+                                        indices[2] = 5;
+                                        break;
+                                case 1:
+                                        indices[0] = 2;
+                                        indices[1] = 0;
+                                        indices[2] = 3;
+                                        break;
+                                case 2:
+                                        indices[0] = 0;
+                                        indices[1] = 1;
+                                        indices[2] = 4;
+                                        break;
+                                default:
+                                        _error_("Edge index provided ("<<index<<") is not between 0 and 2");
+                        }
+                        break;
+                default:
+                        _error_("Element type "<<EnumToStringx(this->element_type)<<" not supported yet");
+        }
+
+        /*Assign output pointer*/
+        *pnumindices = numindices;
+        *pindices    = indices;
+}
+/*}}}*/

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

@@ -55 +55 @@
                 void GetInputDerivativeValue(IssmDouble* pp, IssmDouble* plist,IssmDouble* xyz_list, GaussTria* gauss);
 
+                void NodeOnEdgeIndices(int* pnumindices,int** pindices,int index,int finiteelement);
                 int  NumberofNodes(void);
                 int  NumberofNodes(int finiteelement);