Changeset 16700

Timestamp:

11/11/13 09:13:28 (11 years ago)

Author:

seroussi

Message:

NEW: added vertical solution for stressbalance analysis

File:

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

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

@@ -100 +100 @@
 }/*}}}*/
 void StressbalanceVerticalAnalysis::InputUpdateFromSolution(IssmDouble* solution,Element* element){/*{{{*/
-        _error_("not implemented yet");
-}/*}}}*/
+
+        int          numnodes = element->GetNumberOfNodes();
+        int          numdof=numnodes*1;
+
+        int          i;
+        int          approximation;
+        int*         doflist  = NULL;
+        IssmDouble*  xyz_list = NULL;
+        IssmDouble   rho_ice,g;
+
+        IssmDouble*  values    = xNew<IssmDouble>(numdof);
+        IssmDouble*  vx        = xNew<IssmDouble>(numnodes);
+        IssmDouble*  vy        = xNew<IssmDouble>(numnodes);
+        IssmDouble*  vz        = xNew<IssmDouble>(numnodes);
+        IssmDouble*  vzSSA     = xNew<IssmDouble>(numnodes);
+        IssmDouble*  vzHO      = xNew<IssmDouble>(numnodes);
+        IssmDouble*  vzFS      = xNew<IssmDouble>(numnodes);
+        IssmDouble*  vel       = xNew<IssmDouble>(numnodes);
+        IssmDouble*  pressure  = xNew<IssmDouble>(numnodes);
+        IssmDouble*  surface   = xNew<IssmDouble>(numnodes);
+
+        /*Get the approximation and do nothing if the element in FS or None*/
+        element->GetInputValue(&approximation,ApproximationEnum);
+        if(approximation==FSApproximationEnum || approximation==NoneApproximationEnum){
+                return;
+        }
+
+        /*Get dof list and vertices coordinates: */
+        element->GetVerticesCoordinates(&xyz_list);
+        element->GetDofList(&doflist,NoneApproximationEnum,GsetEnum);
+
+        /*Use the dof list to index into the solution vector vz: */
+        for(i=0;i<numdof;i++) values[i]=solution[doflist[i]];
+        for(i=0;i<numdof;i++){
+                vz[i]=values[i*1+0];
+
+                /*Check solution*/
+                if(xIsNan<IssmDouble>(vz[i])) _error_("NaN found in solution vector");
+        }
+
+        /*Get Vx and Vy*/
+        element->GetInputListOnNodes(&vx[0],VxEnum,0.0); //default is 0
+        element->GetInputListOnNodes(&vy[0],VyEnum,0.0); //default is 0
+
+        /*Do some modifications if we actually have a HOFS or SSAFS element*/
+        if(approximation==HOFSApproximationEnum){
+                Input* vzFS_input=element->GetInput(VzFSEnum);
+                if (vzFS_input){
+                        if (vzFS_input->ObjectEnum()!=PentaInputEnum) _error_("Cannot compute Vel as VzFS is of type " << EnumToStringx(vzFS_input->ObjectEnum()));
+                        element->GetInputListOnNodes(&vzFS[0],VzFSEnum,0.);
+                }
+                else _error_("Cannot compute Vz as VzFS in not present in HOFS element");
+                for(i=0;i<numnodes;i++){
+                        vzHO[i]=vz[i];
+                        vz[i]=vzHO[i]+vzFS[i];
+                }
+        }
+        else if(approximation==SSAFSApproximationEnum){
+                Input* vzFS_input=element->GetInput(VzFSEnum);
+                if (vzFS_input){
+                        if (vzFS_input->ObjectEnum()!=PentaInputEnum) _error_("Cannot compute Vel as VzFS is of type " << EnumToStringx(vzFS_input->ObjectEnum()));
+                        element->GetInputListOnNodes(&vzFS[0],VzFSEnum,0.);
+                }
+                else _error_("Cannot compute Vz as VzFS in not present in SSAFS element");
+                for(i=0;i<numnodes;i++){
+                        vzSSA[i]=vz[i];
+                        vz[i]=vzSSA[i]+vzFS[i];
+                }
+        }
+
+        /*Now Compute vel*/
+        for(i=0;i<numnodes;i++) vel[i]=pow( pow(vx[i],2.0) + pow(vy[i],2.0) + pow(vz[i],2.0) , 0.5);
+
+        /*For pressure: we have not computed pressure in this analysis, for this element. We are in 3D, 
+         *so the pressure is just the pressure at the z elevation: except it this is a HOFS element */
+        if(approximation!=HOFSApproximationEnum &&  approximation!=SSAFSApproximationEnum){
+                rho_ice = element->GetMaterialParameter(MaterialsRhoIceEnum);
+                g       = element->GetMaterialParameter(ConstantsGEnum);
+                element->GetInputListOnNodes(&surface[0],SurfaceEnum,0.);
+                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 Vz inputs to old 
+         * status, otherwise, we'll wipe them off and add the new inputs: */
+        element->InputChangeName(VzEnum,VzPicardEnum);
+
+        if(approximation!=HOFSApproximationEnum && approximation!=SSAFSApproximationEnum){
+                element->InputChangeName(PressureEnum,PressurePicardEnum);
+                element->AddInput(PressureEnum,pressure,P1Enum);
+        }
+        else if(approximation==HOFSApproximationEnum){
+                element->AddInput(VzHOEnum,vzHO,P1Enum);
+        }
+        else if(approximation==SSAFSApproximationEnum){
+                element->AddInput(VzSSAEnum,vzSSA,P1Enum);
+        }
+        element->AddInput(VzEnum,vz,P1Enum);
+        element->AddInput(VelEnum,vel,P1Enum);
+
+        /*Free ressources:*/
+        xDelete<IssmDouble>(surface);
+        xDelete<IssmDouble>(pressure);
+        xDelete<IssmDouble>(vel);
+        xDelete<IssmDouble>(vz);
+        xDelete<IssmDouble>(vzSSA);
+        xDelete<IssmDouble>(vzHO);
+        xDelete<IssmDouble>(vzFS);
+        xDelete<IssmDouble>(vy);
+        xDelete<IssmDouble>(vx);
+        xDelete<IssmDouble>(values);
+        xDelete<IssmDouble>(xyz_list);
+        xDelete<int>(doflist);
+}/*}}}*/