Index: /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp (revision 26203) +++ /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp (revision 26204) @@ -783,7 +783,9 @@ iomodel->FetchDataToInput(inputs,elements,"md.initialization.vx",VxEnum,0.); iomodel->FetchDataToInput(inputs,elements,"md.initialization.vy",VyEnum,0.); - if(isMLHO){//itapopo FIXME the shear velocities should be initialized correctly + if(isMLHO){//itapopo FIXME shear and base velocities should be initialized correctly iomodel->FetchDataToInput(inputs,elements,"md.initialization.vx",VxShearEnum,0.); iomodel->FetchDataToInput(inputs,elements,"md.initialization.vy",VyShearEnum,0.); + iomodel->FetchDataToInput(inputs,elements,"md.initialization.vx",VxBaseEnum,0.); + iomodel->FetchDataToInput(inputs,elements,"md.initialization.vy",VyBaseEnum,0.); } iomodel->FetchDataToInput(inputs,elements,"md.stressbalance.loadingforcex",LoadingforceXEnum); @@ -3202,5 +3204,5 @@ /*Fetch number of nodes and dof for this finite element*/ int numnodes = basalelement->GetNumberOfNodes(); - int numdof = numnodes*2*2; //4 DOFs per node + int numdof = numnodes*4; //4 DOFs per node /*Fetch dof list and allocate solution vectors*/ @@ -3208,8 +3210,10 @@ //basalelement->GetDofListLocal(&doflist,NoneApproximationEnum,GsetEnum); IssmDouble* values = xNew(numdof); - IssmDouble* vx = xNew(numnodes); - IssmDouble* vy = xNew(numnodes); + IssmDouble* vbx = xNew(numnodes); + IssmDouble* vby = xNew(numnodes); IssmDouble* vshx = xNew(numnodes); IssmDouble* vshy = xNew(numnodes); + IssmDouble* vsx = xNew(numnodes); + IssmDouble* vsy = xNew(numnodes); IssmDouble* vz = xNew(numnodes); IssmDouble* vel = xNew(numnodes); @@ -3229,14 +3233,16 @@ /*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(xIsInf(vx[i])) _error_("Inf found in solution vector"); + vsx[i] =vbx[i]+vshx[i]; //surface vx + if(xIsNan(vbx[i])) _error_("NaN found in solution vector"); + if(xIsInf(vbx[i])) _error_("Inf found in solution vector"); if(xIsNan(vshx[i])) _error_("NaN found in solution vector"); if(xIsInf(vshx[i])) _error_("Inf found in solution vector"); - vy[i] =values[i*4+2]; //base vy + vby[i] =values[i*4+2]; //base vy vshy[i]=values[i*4+3]; //shear vy - if(xIsNan(vy[i])) _error_("NaN found in solution vector"); - if(xIsInf(vy[i])) _error_("Inf found in solution vector"); + vsy[i] =vby[i]+vshy[i]; //surface vy + if(xIsNan(vby[i])) _error_("NaN found in solution vector"); + if(xIsInf(vby[i])) _error_("Inf found in solution vector"); if(xIsNan(vshy[i])) _error_("NaN found in solution vector"); if(xIsInf(vshy[i])) _error_("Inf found in solution vector"); @@ -3249,16 +3255,22 @@ /*Get Vz and compute vel (base)*/ basalelement->GetInputListOnNodes(&vz[0],VzEnum,0.); - for(i=0;iAddBasalInput(VxEnum,vx,element->GetElementType()); - element->AddBasalInput(VyEnum,vy,element->GetElementType()); - element->AddBasalInput(VelEnum,vel,element->GetElementType()); + element->AddBasalInput(VxBaseEnum,vbx,element->GetElementType()); + element->AddBasalInput(VyBaseEnum,vby,element->GetElementType()); + element->AddBasalInput(VelEnum,vel,element->GetElementType()); //itapopo check this + + /*Add vx and vy as inputs to the tria element (surface velocities): */ + element->AddBasalInput(VxSurfaceEnum,vsx,element->GetElementType()); + element->AddBasalInput(VySurfaceEnum,vsy,element->GetElementType()); /*Free ressources:*/ xDelete(vel); xDelete(vz); - xDelete(vy); - xDelete(vx); + xDelete(vby); + xDelete(vbx); + xDelete(vsy); + xDelete(vsx); xDelete(vshy); xDelete(vshx); @@ -3270,5 +3282,5 @@ void StressbalanceAnalysis::GetSolutionFromInputsMLHO(Vector* solution,Element* element){/*{{{*/ - IssmDouble vx,vy,vshx,vshy; + IssmDouble vbx,vby,vshx,vshy; int domaintype,dim,approximation,dofpernode; int* doflist = NULL; @@ -3293,7 +3305,7 @@ /*Get inputs*/ - Input* vx_input =element->GetInput(VxEnum); _assert_(vx_input); + Input* vxbase_input =element->GetInput(VxBaseEnum); _assert_(vxbase_input); Input* vxshear_input =element->GetInput(VxShearEnum); _assert_(vxshear_input); - Input* vy_input =element->GetInput(VyEnum); _assert_(vy_input); + Input* vybase_input =element->GetInput(VyBaseEnum); _assert_(vybase_input); Input* vyshear_input =element->GetInput(VyShearEnum); _assert_(vyshear_input); @@ -3304,11 +3316,11 @@ /*Recover vx and vy*/ - vx_input->GetInputValue(&vx,gauss); //base vx + vxbase_input->GetInputValue(&vbx,gauss); //base vx vxshear_input->GetInputValue(&vshx,gauss);//shear vx - values[i*4+0]=vx; //base vx + values[i*4+0]=vbx; //base vx values[i*4+1]=vshx; //shear vx - vy_input->GetInputValue(&vy,gauss); //base vy + vybase_input->GetInputValue(&vby,gauss); //base vy vyshear_input->GetInputValue(&vshy,gauss);//shear vy - values[i*4+2]=vy; //base vy + values[i*4+2]=vby; //base vy values[i*4+3]=vshy;//shear vy } Index: /issm/trunk-jpl/src/c/classes/Loads/Friction.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Loads/Friction.cpp (revision 26203) +++ /issm/trunk-jpl/src/c/classes/Loads/Friction.cpp (revision 26204) @@ -852,9 +852,12 @@ }/*}}}*/ IssmDouble Friction::VelMag(Gauss* gauss){/*{{{*/ - /*Get effective pressure as a function of flag */ - + /*Get the velocity magnitude as a function of flag */ /*diverse*/ + bool isMLHO=false; //employing mono layer higher order model? Default is no IssmDouble vx,vy,vz,vmag; + + element->parameters->FindParam(&isMLHO,FlowequationIsMLHOEnum); + if(isMLHO) _assert_(dim==2); switch(dim){ @@ -864,6 +867,12 @@ break; case 2: - element->GetInputValue(&vx,gauss,VxEnum); - element->GetInputValue(&vy,gauss,VyEnum); + if(isMLHO){ + element->GetInputValue(&vx,gauss,VxBaseEnum); + element->GetInputValue(&vy,gauss,VyBaseEnum); + } + else{ + element->GetInputValue(&vx,gauss,VxEnum); + element->GetInputValue(&vy,gauss,VyEnum); + } vmag=sqrt(vx*vx+vy*vy); break; Index: /issm/trunk-jpl/src/c/shared/Enum/Enum.vim =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/Enum.vim (revision 26203) +++ /issm/trunk-jpl/src/c/shared/Enum/Enum.vim (revision 26204) @@ -936,13 +936,17 @@ syn keyword cConstant VelEnum syn keyword cConstant VxAverageEnum +syn keyword cConstant VxBaseEnum syn keyword cConstant VxEnum syn keyword cConstant VxMeshEnum syn keyword cConstant VxObsEnum syn keyword cConstant VxShearEnum +syn keyword cConstant VxSurfaceEnum syn keyword cConstant VyAverageEnum +syn keyword cConstant VyBaseEnum syn keyword cConstant VyEnum syn keyword cConstant VyMeshEnum syn keyword cConstant VyObsEnum syn keyword cConstant VyShearEnum +syn keyword cConstant VySurfaceEnum syn keyword cConstant VzEnum syn keyword cConstant VzFSEnum Index: /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h (revision 26203) +++ /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h (revision 26204) @@ -933,13 +933,17 @@ VelEnum, VxAverageEnum, + VxBaseEnum, VxEnum, VxMeshEnum, VxObsEnum, VxShearEnum, + VxSurfaceEnum, VyAverageEnum, + VyBaseEnum, VyEnum, VyMeshEnum, VyObsEnum, VyShearEnum, + VySurfaceEnum, VzEnum, VzFSEnum, Index: /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp (revision 26203) +++ /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp (revision 26204) @@ -938,13 +938,17 @@ case VelEnum : return "Vel"; case VxAverageEnum : return "VxAverage"; + case VxBaseEnum : return "VxBase"; case VxEnum : return "Vx"; case VxMeshEnum : return "VxMesh"; case VxObsEnum : return "VxObs"; case VxShearEnum : return "VxShear"; + case VxSurfaceEnum : return "VxSurface"; case VyAverageEnum : return "VyAverage"; + case VyBaseEnum : return "VyBase"; case VyEnum : return "Vy"; case VyMeshEnum : return "VyMesh"; case VyObsEnum : return "VyObs"; case VyShearEnum : return "VyShear"; + case VySurfaceEnum : return "VySurface"; case VzEnum : return "Vz"; case VzFSEnum : return "VzFS"; Index: /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp (revision 26203) +++ /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp (revision 26204) @@ -959,13 +959,17 @@ else if (strcmp(name,"Vel")==0) return VelEnum; else if (strcmp(name,"VxAverage")==0) return VxAverageEnum; + else if (strcmp(name,"VxBase")==0) return VxBaseEnum; else if (strcmp(name,"Vx")==0) return VxEnum; else if (strcmp(name,"VxMesh")==0) return VxMeshEnum; else if (strcmp(name,"VxObs")==0) return VxObsEnum; else if (strcmp(name,"VxShear")==0) return VxShearEnum; + else if (strcmp(name,"VxSurface")==0) return VxSurfaceEnum; else if (strcmp(name,"VyAverage")==0) return VyAverageEnum; + else if (strcmp(name,"VyBase")==0) return VyBaseEnum; else if (strcmp(name,"Vy")==0) return VyEnum; else if (strcmp(name,"VyMesh")==0) return VyMeshEnum; else if (strcmp(name,"VyObs")==0) return VyObsEnum; else if (strcmp(name,"VyShear")==0) return VyShearEnum; + else if (strcmp(name,"VySurface")==0) return VySurfaceEnum; else if (strcmp(name,"Vz")==0) return VzEnum; else if (strcmp(name,"VzFS")==0) return VzFSEnum; @@ -994,12 +998,12 @@ else if (strcmp(name,"Outputdefinition18")==0) return Outputdefinition18Enum; else if (strcmp(name,"Outputdefinition19")==0) return Outputdefinition19Enum; - else if (strcmp(name,"Outputdefinition20")==0) return Outputdefinition20Enum; + else stage=9; + } + if(stage==9){ + if (strcmp(name,"Outputdefinition20")==0) return Outputdefinition20Enum; else if (strcmp(name,"Outputdefinition21")==0) return Outputdefinition21Enum; else if (strcmp(name,"Outputdefinition22")==0) return Outputdefinition22Enum; else if (strcmp(name,"Outputdefinition23")==0) return Outputdefinition23Enum; - else stage=9; - } - if(stage==9){ - if (strcmp(name,"Outputdefinition24")==0) return Outputdefinition24Enum; + else if (strcmp(name,"Outputdefinition24")==0) return Outputdefinition24Enum; else if (strcmp(name,"Outputdefinition25")==0) return Outputdefinition25Enum; else if (strcmp(name,"Outputdefinition26")==0) return Outputdefinition26Enum; @@ -1117,12 +1121,12 @@ else if (strcmp(name,"BoolParam")==0) return BoolParamEnum; else if (strcmp(name,"Boundary")==0) return BoundaryEnum; - else if (strcmp(name,"BuddJacka")==0) return BuddJackaEnum; + else stage=10; + } + if(stage==10){ + if (strcmp(name,"BuddJacka")==0) return BuddJackaEnum; else if (strcmp(name,"CalvingDev2")==0) return CalvingDev2Enum; else if (strcmp(name,"CalvingHab")==0) return CalvingHabEnum; else if (strcmp(name,"CalvingLevermann")==0) return CalvingLevermannEnum; - else stage=10; - } - if(stage==10){ - if (strcmp(name,"CalvingVonmises")==0) return CalvingVonmisesEnum; + else if (strcmp(name,"CalvingVonmises")==0) return CalvingVonmisesEnum; else if (strcmp(name,"Cfdragcoeffabsgrad")==0) return CfdragcoeffabsgradEnum; else if (strcmp(name,"Cfsurfacelogvel")==0) return CfsurfacelogvelEnum; @@ -1240,12 +1244,12 @@ else if (strcmp(name,"IntExternalResult")==0) return IntExternalResultEnum; else if (strcmp(name,"ElementInput")==0) return ElementInputEnum; - else if (strcmp(name,"IntMatExternalResult")==0) return IntMatExternalResultEnum; + else stage=11; + } + if(stage==11){ + if (strcmp(name,"IntMatExternalResult")==0) return IntMatExternalResultEnum; else if (strcmp(name,"IntMatParam")==0) return IntMatParamEnum; else if (strcmp(name,"IntParam")==0) return IntParamEnum; else if (strcmp(name,"IntVecParam")==0) return IntVecParamEnum; - else stage=11; - } - if(stage==11){ - if (strcmp(name,"Inputs")==0) return InputsEnum; + else if (strcmp(name,"Inputs")==0) return InputsEnum; else if (strcmp(name,"Internal")==0) return InternalEnum; else if (strcmp(name,"Intersect")==0) return IntersectEnum; @@ -1363,12 +1367,12 @@ else if (strcmp(name,"SamplingAnalysis")==0) return SamplingAnalysisEnum; else if (strcmp(name,"SamplingSolution")==0) return SamplingSolutionEnum; - else if (strcmp(name,"SIAApproximation")==0) return SIAApproximationEnum; + else stage=12; + } + if(stage==12){ + if (strcmp(name,"SIAApproximation")==0) return SIAApproximationEnum; else if (strcmp(name,"SMBcomponents")==0) return SMBcomponentsEnum; else if (strcmp(name,"SMBd18opdd")==0) return SMBd18opddEnum; else if (strcmp(name,"SMBforcing")==0) return SMBforcingEnum; - else stage=12; - } - if(stage==12){ - if (strcmp(name,"SMBgcm")==0) return SMBgcmEnum; + else if (strcmp(name,"SMBgcm")==0) return SMBgcmEnum; else if (strcmp(name,"SMBgemb")==0) return SMBgembEnum; else if (strcmp(name,"SMBgradients")==0) return SMBgradientsEnum; Index: /issm/trunk-jpl/src/m/solve/parseresultsfromdisk.m =================================================================== --- /issm/trunk-jpl/src/m/solve/parseresultsfromdisk.m (revision 26203) +++ /issm/trunk-jpl/src/m/solve/parseresultsfromdisk.m (revision 26204) @@ -228,4 +228,12 @@ field = field*yts; elseif strcmp(fieldname,'VyShear'), + field = field*yts; + elseif strcmp(fieldname,'VxBase'), + field = field*yts; + elseif strcmp(fieldname,'VyBase'), + field = field*yts; + elseif strcmp(fieldname,'VxSurface'), + field = field*yts; + elseif strcmp(fieldname,'VySurface'), field = field*yts; elseif strcmp(fieldname,'BasalforcingsGroundediceMeltingRate'), Index: /issm/trunk-jpl/src/m/solve/parseresultsfromdisk.py =================================================================== --- /issm/trunk-jpl/src/m/solve/parseresultsfromdisk.py (revision 26203) +++ /issm/trunk-jpl/src/m/solve/parseresultsfromdisk.py (revision 26204) @@ -179,4 +179,12 @@ field = field * yts elif fieldname == 'VyShear': + field = field * yts + elif fieldname == 'VxBase': + field = field * yts + elif fieldname == 'VyBase': + field = field * yts + elif fieldname == 'VxSurface': + field = field * yts + elif fieldname == 'VySurface': field = field * yts elif fieldname == 'BasalforcingsGroundediceMeltingRate':