Changeset 16774

Timestamp:

11/14/13 21:14:19 (11 years ago)

Author:

seroussi

Message:

NEW: working on HOFS

Location:

issm/trunk-jpl/src/c/analyses

Files:

• StressbalanceAnalysis.cpp (modified) (3 diffs)
• StressbalanceAnalysis.h (modified) (1 diff)

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

@@ -950 +950 @@
                         return;
                 case HOFSApproximationEnum:
-                        _error_("not implemented yet");
+                        InputUpdateFromSolutionHOFS(solution,element);
+                        return;
                 case SSAFSApproximationEnum:
                         _error_("not implemented yet");
@@ -1006 +1007 @@
 
         /*Transform solution in Cartesian Space*/
-        element->TransformSolutionCoord(&values[0],&cs_list[0]);
+        element->TransformSolutionCoord(values,cs_list);
 
         /*Ok, we have vx and vy in values, fill in all arrays: */
@@ -1127 +1128 @@
         xDelete<IssmDouble>(xyz_list);
         xDelete<int>(doflist);
+}/*}}}*/
+void StressbalanceAnalysis::InputUpdateFromSolutionHOFS(IssmDouble* solution,Element* element){/*{{{*/
+
+        int         i;
+        IssmDouble  rho_ice,g,FSreconditioning;
+        int*        doflistHO  = NULL;
+        int*        doflistFSv = NULL;
+        int*        doflistFSp = NULL;
+
+        /*Only works with Penta for now*/
+        if(element->ObjectEnum()!=PentaEnum) _error_("Coupling not supported for "<<EnumToStringx(element->ObjectEnum()));
+
+        /*Fetch number of nodes and dof for this finite element*/
+        int numnodes  = 6;
+        int numdofHO  = 6*2;
+        int numdofFSv = 6*3;
+        int numdofFSp = 6;
+
+        /*Fetch dof list and allocate solution vectors*/
+        element->GetDofList(&doflistFSv,FSvelocityEnum,GsetEnum);
+        element->GetDofList(&doflistHO, HOApproximationEnum, GsetEnum);
+        element->GetDofListPressure(&doflistFSp,GsetEnum);
+        IssmDouble* HOvalues  = xNew<IssmDouble>(numdofHO);
+        IssmDouble* FSvalues  = xNew<IssmDouble>(numdofFSv+numdofFSp);
+        IssmDouble* vx        = xNew<IssmDouble>(numnodes);
+        IssmDouble* vy        = xNew<IssmDouble>(numnodes);
+        IssmDouble* vz        = xNew<IssmDouble>(numnodes);
+        IssmDouble* vzHO      = xNew<IssmDouble>(numnodes);
+        IssmDouble* vzFS      = xNew<IssmDouble>(numnodes);
+        IssmDouble* vel       = xNew<IssmDouble>(numnodes);
+        IssmDouble* pressure  = xNew<IssmDouble>(numnodes);
+
+        /*Prepare coordinate system list*/
+        int* cs_list = xNew<int>(2*numnodes);
+        for(i=0;i<numnodes;i++) cs_list[i]          = XYZEnum;
+        for(i=0;i<numnodes;i++) cs_list[numnodes+i] = PressureEnum;
+
+        /*Use the dof list to index into the solution vector: */
+        element->FindParam(&FSreconditioning,StressbalanceFSreconditioningEnum);
+        for(i=0;i<numdofHO ;i++) HOvalues[i]=solution[doflistHO[i]];
+        for(i=0;i<numdofFSv;i++) FSvalues[i]=solution[doflistFSv[i]];
+        for(i=0;i<numdofFSp;i++) FSvalues[numdofFSv+i]=solution[doflistFSp[i]];
+
+        /*Transform solution in Cartesian Space*/
+        element->TransformSolutionCoord(FSvalues,cs_list);
+        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]       = FSvalues[i*3+0]+HOvalues[i*2+0];
+                vy[i]       = FSvalues[i*3+1]+HOvalues[i*2+1];
+                vzFS[i]     = FSvalues[i*3+2];
+                pressure[i] = FSvalues[numnodes*3+i]*FSreconditioning;
+
+                /*Check solution*/
+                if(xIsNan<IssmDouble>(vx[i]))       _error_("NaN found in solution vector");
+                if(xIsNan<IssmDouble>(vy[i]))       _error_("NaN found in solution vector");
+                if(xIsNan<IssmDouble>(vzFS[i]))     _error_("NaN found in solution vector");
+                if(xIsNan<IssmDouble>(pressure[i])) _error_("NaN found in solution vector");
+        }
+
+        /*Get Vz and compute vel*/
+        element->GetInputListOnVertices(vzHO,VzHOEnum);
+        for(i=0;i<numnodes;i++){
+                vz[i] = vzHO[i]+vzFS[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);
+        element->InputChangeName(VzEnum,VzPicardEnum);
+        element->InputChangeName(PressureEnum,PressurePicardEnum);
+
+        /*Add vx and vy as inputs to element: */
+        element->AddInput(VxEnum,vx,P1Enum);
+        element->AddInput(VyEnum,vy,P1Enum);
+        element->AddInput(VzEnum,vz,P1Enum);
+        element->AddInput(VzFSEnum,vzFS,P1Enum);
+        element->AddInput(VelEnum,vel,P1Enum);
+        element->AddInput(PressureEnum,pressure,P1Enum);
+
+        /*Free ressources:*/
+        xDelete<IssmDouble>(pressure);
+        xDelete<IssmDouble>(vel);
+        xDelete<IssmDouble>(vz);
+        xDelete<IssmDouble>(vzHO);
+        xDelete<IssmDouble>(vzFS);
+        xDelete<IssmDouble>(vy);
+        xDelete<IssmDouble>(vx);
+        xDelete<IssmDouble>(FSvalues);
+        xDelete<IssmDouble>(HOvalues);
+        xDelete<int>(doflistFSp);
+        xDelete<int>(doflistFSv);
+        xDelete<int>(doflistHO);
+        xDelete<int>(cs_list);
 }/*}}}*/
 void StressbalanceAnalysis::InputUpdateFromSolutionL1L2(IssmDouble* solution,Element* element){/*{{{*/

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

@@ -24 +24 @@
                 void InputUpdateFromSolutionFS(IssmDouble* solution,Element* element);
                 void InputUpdateFromSolutionHO(IssmDouble* solution,Element* element);
+                void InputUpdateFromSolutionHOFS(IssmDouble* solution,Element* element);
                 void InputUpdateFromSolutionL1L2(IssmDouble* solution,Element* element);
                 void InputUpdateFromSolutionSSA(IssmDouble* solution,Element* element);