source: issm/oecreview/Archive/16554-17801/ISSM-16775-16776.diff@ 17802

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

@@ -26 +26 @@
                 void InputUpdateFromSolutionHOFS(IssmDouble* solution,Element* element);
                 void InputUpdateFromSolutionL1L2(IssmDouble* solution,Element* element);
                 void InputUpdateFromSolutionSSA(IssmDouble* solution,Element* element);
+                void InputUpdateFromSolutionSSAFS(IssmDouble* solution,Element* element);
                 void InputUpdateFromSolutionSSAHO(IssmDouble* solution,Element* element);
 };
 #endif

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

@@ -952 +952 @@
                         InputUpdateFromSolutionHOFS(solution,element);
                         return;
                 case SSAFSApproximationEnum:
-                        _error_("not implemented yet");
+                        InputUpdateFromSolutionSSAFS(solution,element);
                         return;
                 default:
                         _error_("Approximation "<<EnumToStringx(approximation)<<" not supported");
@@ -1423 +1423 @@
         xDelete<int>(doflist);
         if(meshtype!=Mesh2DhorizontalEnum) delete basalelement;
 }/*}}}*/
+void StressbalanceAnalysis::InputUpdateFromSolutionSSAFS(IssmDouble* solution,Element* element){/*{{{*/
+
+        int         i;
+        IssmDouble  rho_ice,g,FSreconditioning;
+        int*        doflistSSA  = NULL;
+        int*        doflistFSv = NULL;
+        int*        doflistFSp = NULL;
+
+        /*we have to add results of this element for FS 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  = 6;
+        int numdof2d  = numnodes;
+        int numdofSSA = 6*2;
+        int numdofFSv = 6*3;
+        int numdofFSp = 6;
+
+        /*Fetch dof list and allocate solution vectors*/
+        element->GetDofList(&doflistFSv,FSvelocityEnum,GsetEnum);
+        element->GetDofListPressure(&doflistFSp,GsetEnum);
+        basalelement->GetDofList(&doflistSSA, SSAApproximationEnum, GsetEnum);
+        IssmDouble* SSAvalues  = xNew<IssmDouble>(numdofSSA);
+        IssmDouble* FSvalues  = xNew<IssmDouble>(numdofFSv+numdofFSp);
+        IssmDouble* vx        = xNew<IssmDouble>(numnodes);
+        IssmDouble* vy        = xNew<IssmDouble>(numnodes);
+        IssmDouble* vz        = xNew<IssmDouble>(numnodes);
+        IssmDouble* vzSSA      = 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<numdof2d;i++){
+                SSAvalues[i]          = solution[doflistSSA[i]];
+                SSAvalues[i+numdof2d] = solution[doflistSSA[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,2*numnodes,cs_list);
+        element->TransformSolutionCoord(SSAvalues,numnodes,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]+SSAvalues[i*2+0];
+                vy[i]       = FSvalues[i*3+1]+SSAvalues[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(vzSSA,VzSSAEnum);
+        for(i=0;i<numnodes;i++){
+                vz[i] = vzSSA[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>(vzSSA);
+        xDelete<IssmDouble>(vzFS);
+        xDelete<IssmDouble>(vy);
+        xDelete<IssmDouble>(vx);
+        xDelete<IssmDouble>(FSvalues);
+        xDelete<IssmDouble>(SSAvalues);
+        xDelete<int>(doflistFSp);
+        xDelete<int>(doflistFSv);
+        xDelete<int>(doflistSSA);
+        xDelete<int>(cs_list);
+}/*}}}*/
 void StressbalanceAnalysis::InputUpdateFromSolutionSSAHO(IssmDouble* solution,Element* element){/*{{{*/
 
         int         i,meshtype;