Changeset 16727

Timestamp:

11/12/13 20:38:11 (11 years ago)

Author:

seroussi

Message:

BUG: fixed SSA InputUpdateFromSolution

Location:

issm/trunk-jpl/src/c

Files:

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

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

@@ -297 +297 @@
         xDelete<IssmDouble>(phi);
         xDelete<int>(doflist);
-        if(meshtype!=Mesh2DhorizontalEnum) delete element;
-}/*}}}*/
+        if(meshtype!=Mesh2DhorizontalEnum) delete basalelement;
+}/*}}}*/

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

@@ -937 +937 @@
                         InputUpdateFromSolutionFS(solution,element);
                         return;
-                case SSAApproximationEnum: case HOApproximationEnum: case SIAApproximationEnum:
-                        InputUpdateFromSolutionHoriz(solution,element);
+                case SSAApproximationEnum: 
+                        InputUpdateFromSolutionSSA(solution,element);
+                        return;
+                case HOApproximationEnum: 
+                        InputUpdateFromSolutionHO(solution,element);
+                        return;
+                case SIAApproximationEnum:
+                        InputUpdateFromSolutionHO(solution,element);
                         return;
                 case L1L2ApproximationEnum:
-                        InputUpdateFromSolutionHoriz(solution,element);
+                        InputUpdateFromSolutionHO(solution,element);
                         return;
                 case SSAHOApproximationEnum: case HOFSApproximationEnum: case SSAFSApproximationEnum:
@@ -950 +956 @@
         }
 }/*}}}*/
-void StressbalanceAnalysis::InputUpdateFromSolutionHoriz(IssmDouble* solution,Element* element){/*{{{*/
+void StressbalanceAnalysis::InputUpdateFromSolutionSSA(IssmDouble* solution,Element* element){/*{{{*/
+
+        int         i,meshtype;
+        IssmDouble  rho_ice,g;
+        int*        doflist=NULL;
+        IssmDouble* xyz_list=NULL;
+        Element*    basalelement=NULL;
+
+        /*Deal with pressure first*/
+        int numvertices = element->GetNumberOfVertices();
+        IssmDouble* pressure  = xNew<IssmDouble>(numvertices);
+        IssmDouble* thickness = xNew<IssmDouble>(numvertices);
+        IssmDouble* surface   = xNew<IssmDouble>(numvertices);
+
+        element->FindParam(&meshtype,MeshTypeEnum);
+        rho_ice =element->GetMaterialParameter(MaterialsRhoIceEnum);
+        g       =element->GetMaterialParameter(ConstantsGEnum);
+        switch(meshtype){
+                case Mesh2DhorizontalEnum:
+                        element->GetInputListOnNodes(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);
+                        for(i=0;i<numvertices;i++) pressure[i]=rho_ice*g*(surface[i]-xyz_list[i*3+2]);
+                        break;
+                default: _error_("mesh "<<EnumToStringx(meshtype)<<" not supported yet");
+        }
+        element->AddInput(PressureEnum,pressure,P1Enum);
+        xDelete<IssmDouble>(pressure);
+        xDelete<IssmDouble>(thickness);
+        xDelete<IssmDouble>(surface);
+
+        /*Get basal element*/
+        switch(meshtype){
+                case Mesh2DhorizontalEnum:
+                        basalelement = element;
+                        break;
+                case Mesh3DEnum:
+                        if(!element->IsOnBed()) return;
+                        basalelement=element->SpawnBasalElement();
+                        break;
+                default: _error_("mesh "<<EnumToStringx(meshtype)<<" not supported yet");
+        }
+
+        /*Fetch number of nodes and dof for this finite element*/
+        int numnodes = basalelement->GetNumberOfNodes();
+        int numdof   = numnodes*2;
+
+        /*Fetch dof list and allocate solution vectors*/
+        basalelement->GetDofList(&doflist,NoneApproximationEnum,GsetEnum);
+        IssmDouble* values    = xNew<IssmDouble>(numdof);
+        IssmDouble* vx        = xNew<IssmDouble>(numnodes);
+        IssmDouble* vy        = xNew<IssmDouble>(numnodes);
+        IssmDouble* vz        = xNew<IssmDouble>(numnodes);
+        IssmDouble* vel       = xNew<IssmDouble>(numnodes);
+
+        /*Use the dof list to index into the solution vector: */
+        for(i=0;i<numdof;i++) values[i]=solution[doflist[i]];
+
+        /*Transform solution in Cartesian Space*/
+        basalelement->TransformSolutionCoord(&values[0],XYEnum);
+        basalelement->FindParam(&meshtype,MeshTypeEnum);
+
+        /*Ok, we have vx and vy in values, fill in vx and vy arrays: */
+        for(i=0;i<numnodes;i++){
+                vx[i]=values[i*2+0];
+                vy[i]=values[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*/
+        basalelement->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]);
+
+        /*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);
+
+        /*Add vx and vy as inputs to the tria element: */
+        element->AddBasalInput(VxEnum,vx,P1Enum);
+        element->AddBasalInput(VyEnum,vy,P1Enum);
+        element->AddBasalInput(VelEnum,vel,P1Enum);
+
+        /*Free ressources:*/
+        xDelete<IssmDouble>(vel);
+        xDelete<IssmDouble>(vz);
+        xDelete<IssmDouble>(vy);
+        xDelete<IssmDouble>(vx);
+        xDelete<IssmDouble>(values);
+        xDelete<IssmDouble>(xyz_list);
+        xDelete<int>(doflist);
+        if(meshtype!=Mesh2DhorizontalEnum) delete basalelement;
+}/*}}}*/
+void StressbalanceAnalysis::InputUpdateFromSolutionHO(IssmDouble* solution,Element* element){/*{{{*/
 
         int         i,meshtype;
@@ -969 +1074 @@
         IssmDouble* vel       = xNew<IssmDouble>(numnodes);
         IssmDouble* pressure  = xNew<IssmDouble>(numnodes);
-        IssmDouble* thickness = xNew<IssmDouble>(numnodes);
         IssmDouble* surface   = xNew<IssmDouble>(numnodes);
 
@@ -997 +1101 @@
         rho_ice = element->GetMaterialParameter(MaterialsRhoIceEnum);
         g       = element->GetMaterialParameter(ConstantsGEnum);
-        switch(meshtype){
-                case Mesh2DhorizontalEnum:
-                        element->GetInputListOnNodes(&thickness[0],ThicknessEnum);
-                        for(i=0;i<numnodes;i++) pressure[i]=rho_ice*g*thickness[i];
-                        break;
-                case Mesh3DEnum:   
-                        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]);
-                        break;
-                default: _error_("mesh "<<EnumToStringx(meshtype)<<" not supported yet");
-        }
+        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 
@@ -1019 +1114 @@
         element->AddInput(VxEnum,vx,P1Enum);
         element->AddInput(VyEnum,vy,P1Enum);
-        element->AddInput(VelEnum,vel,P1Enum);
-        element->AddInput(PressureEnum,pressure,P1Enum);
+        element->AddBasalInput(VelEnum,vel,P1Enum);
+        element->AddBasalInput(PressureEnum,pressure,P1Enum);
 
         /*Free ressources:*/
-        xDelete<IssmDouble>(thickness);
         xDelete<IssmDouble>(surface);
         xDelete<IssmDouble>(pressure);
@@ -1033 +1127 @@
         xDelete<IssmDouble>(xyz_list);
         xDelete<int>(doflist);
-
 }/*}}}*/
 void StressbalanceAnalysis::InputUpdateFromSolutionFS(IssmDouble* solution,Element* element){/*{{{*/

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

@@ -23 +23 @@
                 void InputUpdateFromSolution(IssmDouble* solution,Element* element);
                 void InputUpdateFromSolutionFS(IssmDouble* solution,Element* element);
-                void InputUpdateFromSolutionHoriz(IssmDouble* solution,Element* element);
+                void InputUpdateFromSolutionHO(IssmDouble* solution,Element* element);
+                void InputUpdateFromSolutionSSA(IssmDouble* solution,Element* element);
 };
 #endif

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

@@ -60 +60 @@
                 virtual int    GetNumberOfNodesVelocity(void)=0;
                 virtual int    GetNumberOfNodesPressure(void)=0;
+                virtual int    GetNumberOfVertices(void)=0;
                 virtual void   GetNodesSidList(int* sidlist)=0;
                 virtual void   GetNodesLidList(int* sidlist)=0;

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

@@ -139 +139 @@
                         Penta*     penta=NULL;
 
-                        for(i=1;i<NUMVERTICES2D;i++){
+                        for(i=0;i<NUMVERTICES2D;i++){
                                 extrudedvalues[i]=values[i];
                                 extrudedvalues[i+NUMVERTICES2D]=values[i];
                         }
-                        this->inputs->AddInput(new PentaInput(input_enum,&extrudedvalues[0],P1Enum));
                         penta=this;
                         for(;;){
@@ -1394 +1393 @@
 int Penta::GetNumberOfNodesVelocity(void){
         return this->NumberofNodesVelocity();
+}
+/*}}}*/
+/*FUNCTION Penta::GetNumberOfVertices;{{{*/
+int Penta::GetNumberOfVertices(void){
+        return NUMVERTICES; 
 }
 /*}}}*/

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

@@ -93 +93 @@
                 int    GetNumberOfNodesPressure(void);
                 int    GetNumberOfNodesVelocity(void);
+                int    GetNumberOfVertices(void);
                 IssmDouble GetMaterialParameter(int enum_in);
                 void   GetSolutionFromInputsOneDof(Vector<IssmDouble>* solution,int enum_type);

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

@@ -98 +98 @@
                 int         GetNumberOfNodesVelocity(void){_error_("not implemented yet");};
                 int         GetNumberOfNodesPressure(void){_error_("not implemented yet");};
+                int         GetNumberOfVertices(void){_error_("not implemented yet");};
                 void        GetVerticesCoordinates(IssmDouble** pxyz_list){_error_("not implemented yet");};
                 int         Sid(){_error_("not implemented yet");};

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

@@ -1367 +1367 @@
 int Tria::GetNumberOfNodes(void){
         return this->NumberofNodes();
+}
+/*}}}*/
+/*FUNCTION Tria::GetNumberOfVertices;{{{*/
+int Tria::GetNumberOfVertices(void){
+        return NUMVERTICES;
 }
 /*}}}*/

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

@@ -91 +91 @@
                 int         GetNumberOfNodesPressure(void){_error_("not implemented yet");};
                 int         GetNumberOfNodesVelocity(void){_error_("not implemented yet");};
+                int         GetNumberOfVertices(void);
                 int         Sid();
                 bool        IsOnBed();