Changeset 16762

Timestamp:

11/14/13 12:15:52 (11 years ago)

Author:

Mathieu Morlighem

Message:

BUG: fixed issues for higher order elements and trying to get Tiling HO SSA to work

Location:

issm/trunk-jpl/src/c

Files:

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

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

@@ -947 +947 @@
                         return;
                 case L1L2ApproximationEnum:
-                        InputUpdateFromSolutionHO(solution,element);
+                        InputUpdateFromSolutionSSA(solution,element);
                         return;
-                case SSAHOApproximationEnum: case HOFSApproximationEnum: case SSAFSApproximationEnum:
-                        /*the elements around will create the solution*/
+                case SSAHOApproximationEnum:
+                        InputUpdateFromSolutionSSAHO(solution,element);
+                        return;
+                case HOFSApproximationEnum:
+                        _error_("not implemented yet");
+                case SSAFSApproximationEnum:
+                        _error_("not implemented yet");
                         return;
                 default:
@@ -975 +980 @@
         switch(meshtype){
                 case Mesh2DhorizontalEnum:
-                        element->GetInputListOnNodes(thickness,ThicknessEnum);
+                        element->GetInputListOnVertices(thickness,ThicknessEnum);
                         for(i=0;i<numvertices;i++) pressure[i]=rho_ice*g*thickness[i];
                         break;
                 case Mesh3DEnum:
                         element->GetVerticesCoordinates(&xyz_list);
-                        element->GetInputListOnNodes(surface,SurfaceEnum);
+                        element->GetInputListOnVertices(surface,SurfaceEnum);
                         for(i=0;i<numvertices;i++) pressure[i]=rho_ice*g*(surface[i]-xyz_list[i*3+2]);
                         break;
@@ -996 +1001 @@
                         break;
                 case Mesh3DEnum:
-                        if(!element->IsOnBed()) return;
+                        if(!element->IsOnBed()){xDelete<IssmDouble>(xyz_list); return;}
                         basalelement=element->SpawnBasalElement();
                         break;
@@ -1225 +1230 @@
         xDelete<int>(cs_list);
 }/*}}}*/
+void StressbalanceAnalysis::InputUpdateFromSolutionSSAHO(IssmDouble* solution,Element* element){/*{{{*/
+
+        int         i,meshtype;
+        IssmDouble  rho_ice,g;
+        int*        SSAdoflist = NULL;
+        int*        HOdoflist  = NULL;
+        IssmDouble* xyz_list   = NULL;
+
+        /*we have to add results of this element for HO and results from the element
+         * at base for SSA, so we need to have the pointer toward the basal element*/
+        Element* basalelement=element->GetBasalElement();
+        if(basalelement->ObjectEnum()!=PentaEnum){
+                _error_("Coupling not supported for "<<EnumToStringx(basalelement->ObjectEnum()));
+        }
+
+        /*Fetch number of nodes and dof for this finite element*/
+        int numnodes = element->GetNumberOfNodes();
+        int numdof   = numnodes*2;
+        int numdof2d = numnodes;/*FIXME: only works for Penta*/
+
+        /*Fetch dof list and allocate solution vectors*/
+        basalelement->GetDofList(&SSAdoflist,SSAApproximationEnum,GsetEnum);
+        element->GetDofList(     &HOdoflist, HOApproximationEnum, GsetEnum);
+        IssmDouble* HOvalues  = xNew<IssmDouble>(numdof);
+        IssmDouble* SSAvalues = xNew<IssmDouble>(numdof);
+        IssmDouble* vx        = xNew<IssmDouble>(numnodes);
+        IssmDouble* vy        = xNew<IssmDouble>(numnodes);
+        IssmDouble* vz        = xNew<IssmDouble>(numnodes);
+        IssmDouble* vel       = xNew<IssmDouble>(numnodes);
+        IssmDouble* pressure  = xNew<IssmDouble>(numnodes);
+        IssmDouble* surface   = xNew<IssmDouble>(numnodes);
+
+        /*Use the dof list to index into the solution vector: */
+        for(i=0;i<numdof2d;i++){
+                HOvalues[i]  = solution[HOdoflist[i]];
+                SSAvalues[i] = solution[SSAdoflist[i]];
+        }
+        for(i=numdof2d;i<numdof;i++){
+                HOvalues[i]  = solution[HOdoflist[i]];
+                SSAvalues[i] = SSAvalues[i-numdof2d];
+        }
+
+        /*Transform solution in Cartesian Space*/
+        basalelement->TransformSolutionCoord(SSAvalues,XYEnum);
+        element->TransformSolutionCoord(HOvalues,XYEnum);
+
+        /*Ok, we have vx and vy in values, fill in vx and vy arrays: */
+        for(i=0;i<numnodes;i++){
+                vx[i]=SSAvalues[i*2+0]+HOvalues[i*2+0];
+                vy[i]=SSAvalues[i*2+1]+HOvalues[i*2+1];
+
+                /*Check solution*/
+                if(xIsNan<IssmDouble>(vx[i])) _error_("NaN found in solution vector");
+                if(xIsNan<IssmDouble>(vy[i])) _error_("NaN found in solution vector");
+        }
+
+        /*Get Vz and compute vel*/
+        element->GetInputListOnNodes(&vz[0],VzEnum,0.);
+        for(i=0;i<numnodes;i++) vel[i]=sqrt(vx[i]*vx[i] + vy[i]*vy[i] + vz[i]*vz[i]);
+
+        /*For pressure: we have not computed pressure in this analysis, for this element. We are in 2D, 
+         *so the pressure is just the pressure at the bedrock: */
+        rho_ice = element->GetMaterialParameter(MaterialsRhoIceEnum);
+        g       = element->GetMaterialParameter(ConstantsGEnum);
+        element->GetVerticesCoordinates(&xyz_list);
+        element->GetInputListOnNodes(&surface[0],SurfaceEnum);
+        for(i=0;i<numnodes;i++) pressure[i]=rho_ice*g*(surface[i]-xyz_list[i*3+2]);
+
+        /*Now, we have to move the previous Vx and Vy inputs  to old 
+         * status, otherwise, we'll wipe them off: */
+        element->InputChangeName(VxEnum,VxPicardEnum);
+        element->InputChangeName(VyEnum,VyPicardEnum);
+        element->InputChangeName(PressureEnum,PressurePicardEnum);
+
+        /*Add vx and vy as inputs to the tria element: */
+        element->AddInput(VxEnum,vx,P1Enum);
+        element->AddInput(VyEnum,vy,P1Enum);
+        element->AddBasalInput(VelEnum,vel,P1Enum);
+        element->AddBasalInput(PressureEnum,pressure,P1Enum);
+
+        /*Free ressources:*/
+        xDelete<IssmDouble>(surface);
+        xDelete<IssmDouble>(pressure);
+        xDelete<IssmDouble>(vel);
+        xDelete<IssmDouble>(vz);
+        xDelete<IssmDouble>(vy);
+        xDelete<IssmDouble>(vx);
+        xDelete<IssmDouble>(xyz_list);
+        xDelete<IssmDouble>(SSAvalues);
+        xDelete<IssmDouble>(HOvalues);
+        xDelete<int>(SSAdoflist);
+        xDelete<int>(HOdoflist);
+}/*}}}*/

issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.h ¶

@@ -25 +25 @@
                 void InputUpdateFromSolutionHO(IssmDouble* solution,Element* element);
                 void InputUpdateFromSolutionSSA(IssmDouble* solution,Element* element);
+                void InputUpdateFromSolutionSSAHO(IssmDouble* solution,Element* element);
 };
 #endif

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

@@ -53 +53 @@
                 virtual void   TransformSolutionCoord(IssmDouble* values,int transformenum)=0;
                 virtual void   TransformSolutionCoord(IssmDouble* values,int* transformenum_list)=0;
+                virtual Element* GetBasalElement(void)=0;
                 virtual void    GetDofList(int** pdoflist,int approximation_enum,int setenum)=0;
                 virtual void    GetDofListVelocity(int** pdoflist,int setenum)=0;

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

@@ -1029 +1029 @@
 /*}}}*/
 /*FUNCTION Penta::GetBasalElement{{{*/
-Penta* Penta::GetBasalElement(void){
+Element* Penta::GetBasalElement(void){
 
         /*Output*/
@@ -5215 +5215 @@
         
         /*get basal element, needed for basal watercolumn*/
-        pentabase=this->GetBasalElement();
+        pentabase=(Penta*)this->GetBasalElement();
         
         this->GetInputListOnVertices(&enthalpy[0],EnthalpyEnum);
@@ -6988 +6988 @@
 
         /*Find penta on bed as HO must be coupled to the dofs on the bed: */
-        Penta* pentabase=GetBasalElement();
+        Penta* pentabase=(Penta*)GetBasalElement();
         Tria*  tria=pentabase->SpawnTria(0); //lower face is 0, upper face is 1.
 
@@ -7183 +7183 @@
 
         /*Find penta on bed as FS must be coupled to the dofs on the bed: */
-        Penta* pentabase=GetBasalElement();
+        Penta* pentabase=(Penta*)GetBasalElement();
         Tria* tria=pentabase->SpawnTria(0); //lower face is 0, upper face is 1.
 
@@ -7618 +7618 @@
 
         /*Find penta on bed as this is a SSA elements: */
-        pentabase=GetBasalElement();
+        pentabase=(Penta*)GetBasalElement();
         tria=pentabase->SpawnTria(0); //lower face is 0, upper face is 1.
 
@@ -7770 +7770 @@
 
         /*Find penta on bed as this is a SSA elements: */
-        pentabase=GetBasalElement();
+        pentabase=(Penta*)GetBasalElement();
         tria=pentabase->SpawnTria(0); //lower face is 0, upper face is 1.
 
@@ -9982 +9982 @@
         /*OK, we have to add results of this element for HO 
          * and results from the penta at base for SSA. Now recover results*/
-        penta=GetBasalElement();
+        penta=(Penta*)GetBasalElement();
 
         /*Get dof listof this element (HO dofs) and of the penta at base (SSA dofs): */
@@ -10075 +10075 @@
         /*OK, we have to add results of this element for SSA 
          * and results from the penta at base for SSA. Now recover results*/
-        penta=GetBasalElement();
+        penta=(Penta*)GetBasalElement();
 
         /*Get dof listof this element (SSA dofs) and of the penta at base (SSA dofs): */

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

@@ -85 +85 @@
                 void   Delta18oParameterization(void);
                 void   EnthalpyToThermal(IssmDouble* ptemperature,IssmDouble* pwaterfraction,IssmDouble enthalpy,IssmDouble pressure);
+                Element* GetBasalElement(void);
                 void     GetDofList(int** pdoflist,int approximation_enum,int setenum);
                 void     GetDofListVelocity(int** pdoflist,int setenum);
@@ -228 +229 @@
                 void           GetZeroLevelsetCoordinates(IssmDouble* xyz_zero,IssmDouble xyz_list[6][3],int levelsetenum);
                 Penta*         GetLowerElement(void);
-                Penta*         GetBasalElement(void);
                 void           InputChangeName(int input_enum, int enum_type_old);
                 void             InputExtrude(int enum_type,int object_type);

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

@@ -87 +87 @@
                 void        FindParam(IssmDouble* pvalue,int paramenum);
                 int         FiniteElement(void);
+                Element*    GetBasalElement(void){_error_("not implemented yet");};
                 void        GetDofList(int** pdoflist,int approximation_enum,int setenum){_error_("not implemented yet");};
                 void        GetDofListVelocity(int** pdoflist,int setenum){_error_("not implemented yet");};

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

@@ -83 +83 @@
                 void        FindParam(IssmDouble* pvalue,int paramenum);
                 int         FiniteElement(void);
+                Element*    GetBasalElement(void){_error_("not implemented yet");};
                 void          GetDofList(int** pdoflist,int approximation_enum,int setenum);
                 void          GetDofListVelocity(int** pdoflist,int setenum);