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 "<(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(numdof); + IssmDouble* vx = xNew(numdof); + IssmDouble* vy = xNew(numdof); + IssmDouble* vz = xNew(numdof); + IssmDouble* vel = xNew(numdof); + IssmDouble* pressure = xNew(numdof); + IssmDouble* thickness = xNew(numdof); + IssmDouble* surface = xNew(numnodes); + + /*Use the dof list to index into the solution vector: */ + for(i=0;iTransformSolutionCoord(&values[0],XYEnum); + + /*Ok, we have vx and vy in values, fill in vx and vy arrays: */ + for(i=0;i(vx[i])) _error_("NaN found in solution vector"); + if(xIsNan(vy[i])) _error_("NaN found in solution vector"); + } + + /*Get Vz and compute vel*/ + element->GetInputListOnNodes(&vz[0],VzEnum,0.); + for(i=0;iGetMaterialParameter(MaterialsRhoIceEnum); + g = element->GetMaterialParameter(ConstantsGEnum); + switch(meshtype){ + case Mesh2DhorizontalEnum: + element->GetInputListOnNodes(&thickness[0],ThicknessEnum); + for(i=0;iGetVerticesCoordinates(&xyz_list); + element->GetInputListOnNodes(&surface[0],SurfaceEnum); + for(i=0;iInputChangeName(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(thickness); + xDelete(surface); + xDelete(pressure); + xDelete(vel); + xDelete(vz); + xDelete(vy); + xDelete(vx); + xDelete(values); + xDelete(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");};