source: issm/oecreview/Archive/16554-17801/ISSM-16716-16717.diff

Index: ../trunk-jpl/src/c/analyses/L2ProjectionBaseAnalysis.cpp
===================================================================
--- ../trunk-jpl/src/c/analyses/L2ProjectionBaseAnalysis.cpp    (revision 16716)
+++ ../trunk-jpl/src/c/analyses/L2ProjectionBaseAnalysis.cpp    (revision 16717)
@@ -59,8 +59,17 @@
 }/*}}}*/
 void L2ProjectionBaseAnalysis::InputUpdateFromSolution(IssmDouble* solution,Element* element){/*{{{*/
 
-       int inputenum;
+       int inputenum,meshtype;
 
        element->FindParam(&inputenum,InputToL2ProjectEnum);
-       element->InputUpdateFromSolutionOneDofCollapsed(solution,inputenum);
+       element->FindParam(&meshtype,MeshTypeEnum);
+       switch(meshtype){
+               case Mesh2DhorizontalEnum:
+                       element->InputUpdateFromSolutionOneDof(solution,inputenum);
+                       break;
+               case Mesh3DEnum:
+                       element->InputUpdateFromSolutionOneDofCollapsed(solution,inputenum);
+                       break;
+               default: _error_("mesh "<<EnumToStringx(meshtype)<<" not supported yet");
+       }
 }/*}}}*/
Index: ../trunk-jpl/src/c/analyses/StressbalanceSIAAnalysis.cpp
===================================================================
--- ../trunk-jpl/src/c/analyses/StressbalanceSIAAnalysis.cpp    (revision 16716)
+++ ../trunk-jpl/src/c/analyses/StressbalanceSIAAnalysis.cpp    (revision 16717)
@@ -167,5 +167,84 @@
        xDelete<IssmDouble>(values);
 }/*}}}*/
 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);
 }/*}}}*/
Index: ../trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp
===================================================================
--- ../trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp       (revision 16716)
+++ ../trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp       (revision 16717)
@@ -940,7 +940,7 @@
                        InputUpdateFromSolutionHoriz(solution,element);
                        return;
                case L1L2ApproximationEnum:
-                       //InputUpdateFromSolutionHoriz(solution,element);
+                       InputUpdateFromSolutionHoriz(solution,element);
                        return;
                case SSAHOApproximationEnum: case HOFSApproximationEnum: case SSAFSApproximationEnum:
                        /*the elements around will create the solution*/
Index: ../trunk-jpl/src/c/classes/Elements/Element.h
===================================================================
--- ../trunk-jpl/src/c/classes/Elements/Element.h       (revision 16716)
+++ ../trunk-jpl/src/c/classes/Elements/Element.h       (revision 16717)
@@ -88,6 +88,7 @@
                virtual void   InputDuplicate(int original_enum,int new_enum)=0;
                virtual void   InputCreate(IssmDouble* vector,IoModel* iomodel,int M,int N,int vector_type,int vector_enum,int code)=0;
                virtual void   InputUpdateFromSolutionOneDofCollapsed(IssmDouble* solution,int inputenum)=0;
+               virtual void   InputUpdateFromSolutionOneDof(IssmDouble* solution,int inputenum)=0;
 
                virtual int    NodalValue(IssmDouble* pvalue, int index, int natureofdataenum)=0;
                virtual int    NumberofNodesVelocity(void)=0;
Index: ../trunk-jpl/src/c/classes/Elements/Seg.h
===================================================================
--- ../trunk-jpl/src/c/classes/Elements/Seg.h   (revision 16716)
+++ ../trunk-jpl/src/c/classes/Elements/Seg.h   (revision 16717)
@@ -55,6 +55,7 @@
                /*Update virtual functions resolution: {{{*/
                void  InputUpdateFromSolution(IssmDouble* solution){_error_("not implemented yet");};
                void  InputUpdateFromSolutionOneDofCollapsed(IssmDouble* solution,int inputenum){_error_("not implemented yet");};
+               void  InputUpdateFromSolutionOneDof(IssmDouble* solution,int inputenum){_error_("not implemented yet");};
                void  InputUpdateFromVector(IssmDouble* vector, int name, int type){_error_("not implemented yet");};
 #ifdef _HAVE_DAKOTA_
                void  InputUpdateFromVectorDakota(IssmDouble* vector, int name, int type){_error_("not implemented yet");};