Changeset 16717 for issm/trunk-jpl/src/c/analyses/StressbalanceSIAAnalysis.cpp

Timestamp:

11/12/13 13:54:57 (11 years ago)

Author:

seroussi

Message:

NEW: moved SIA update from solution in analyses

File:

• issm/trunk-jpl/src/c/analyses/StressbalanceSIAAnalysis.cpp (modified) (1 diff)

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

@@ -168 +168 @@
 }/*}}}*/
 void StressbalanceSIAAnalysis::InputUpdateFromSolution(IssmDouble* solution,Element* element){/*{{{*/
-        _error_("not implemented yet");
-}/*}}}*/
+
+        int         i,meshtype;
+        IssmDouble  rho_ice,g;
+        int*        doflist=NULL;
+        IssmDouble* xyz_list=NULL;
+
+        /*Fetch number of nodes and dof for this finite element*/
+        int numnodes = element->GetNumberOfNodes();
+        int numdof   = numnodes*2;
+
+        /*Fetch dof list and allocate solution vectors*/
+        element->GetDofList(&doflist,NoneApproximationEnum,GsetEnum);
+        IssmDouble* values    = xNew<IssmDouble>(numdof);
+        IssmDouble* vx        = xNew<IssmDouble>(numdof);
+        IssmDouble* vy        = xNew<IssmDouble>(numdof);
+        IssmDouble* vz        = xNew<IssmDouble>(numdof);
+        IssmDouble* vel       = xNew<IssmDouble>(numdof);
+        IssmDouble* pressure  = xNew<IssmDouble>(numdof);
+        IssmDouble* thickness = xNew<IssmDouble>(numdof);
+        IssmDouble* surface   = 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*/
+        element->TransformSolutionCoord(&values[0],XYEnum);
+
+        /*Ok, we have vx and vy in values, fill in vx and vy arrays: */
+        for(i=0;i<numnodes;i++){
+                vx[i]=values[i*NDOF2+0];
+                vy[i]=values[i*NDOF2+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);
+        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");
+        }
+
+        /*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->AddInput(VelEnum,vel,P1Enum);
+        element->AddInput(PressureEnum,pressure,P1Enum);
+
+        /*Free ressources:*/
+        xDelete<IssmDouble>(thickness);
+        xDelete<IssmDouble>(surface);
+        xDelete<IssmDouble>(pressure);
+        xDelete<IssmDouble>(vel);
+        xDelete<IssmDouble>(vz);
+        xDelete<IssmDouble>(vy);
+        xDelete<IssmDouble>(vx);
+        xDelete<IssmDouble>(values);
+        xDelete<int>(doflist);
+}/*}}}*/