Index: /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp (revision 26191) +++ /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp (revision 26192) @@ -194,10 +194,10 @@ } } - else{//itapopo testing here + else{//itapopo FIXME only base velocities are constrained here IoModelToConstraintsx(constraints,iomodel,"md.stressbalance.spcvx",StressbalanceAnalysisEnum,finiteelement,0); - IoModelToConstraintsx(constraints,iomodel,"md.stressbalance.spcvx",StressbalanceAnalysisEnum,finiteelement,1); + //IoModelToConstraintsx(constraints,iomodel,"md.stressbalance.spcvx",StressbalanceAnalysisEnum,finiteelement,1); if(iomodel->domaintype!=Domain2DverticalEnum){ IoModelToConstraintsx(constraints,iomodel,"md.stressbalance.spcvy",StressbalanceAnalysisEnum,finiteelement,2); - IoModelToConstraintsx(constraints,iomodel,"md.stressbalance.spcvy",StressbalanceAnalysisEnum,finiteelement,3); + //IoModelToConstraintsx(constraints,iomodel,"md.stressbalance.spcvy",StressbalanceAnalysisEnum,finiteelement,3); } } @@ -783,7 +783,7 @@ iomodel->FetchDataToInput(inputs,elements,"md.initialization.vx",VxEnum,0.); iomodel->FetchDataToInput(inputs,elements,"md.initialization.vy",VyEnum,0.); - if(isMLHO){//itapopo - iomodel->FetchDataToInput(inputs,elements,"md.initialization.vx",VxBaseEnum,0.); - iomodel->FetchDataToInput(inputs,elements,"md.initialization.vy",VyBaseEnum,0.); + if(isMLHO){//itapopo FIXME the shear 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.stressbalance.loadingforcex",LoadingforceXEnum); @@ -2084,5 +2084,5 @@ /*Use the dof list to index into the solution vector: */ for(i=0;iTransformSolutionCoord(&values[0],XYEnum); @@ -2729,5 +2729,7 @@ /*MLHO*/ ElementMatrix* StressbalanceAnalysis::CreateKMatrixMLHO(Element* element){/*{{{*/ - + + _error_("Mono Layer Higher-Order called, not fully tested. If you are sure in using it, comment this line."); + /* Check if ice in element */ if(!element->IsIceInElement()) return NULL; @@ -2745,7 +2747,87 @@ if(!element->IsOnBase() || element->IsFloating()) return NULL; + + /*Intermediaries*/ + int dim,domaintype; + bool mainlyfloating; + int friction_style,point1; + IssmDouble alpha2,Jdet,fraction1,fraction2; + IssmDouble gllevelset,phi=1.; + IssmDouble *xyz_list = NULL; + Gauss* gauss = NULL; + + /*Get problem dimension*/ + element->FindParam(&domaintype,DomainTypeEnum); + switch(domaintype){ + case Domain2DverticalEnum: dim = 1;break; + case Domain2DhorizontalEnum: dim = 2;break; + case Domain3DEnum: dim = 2;break; + default: _error_("mesh "<GetNumberOfNodes(); + int numdof = numnodes*dim; + + /*Initialize Element matrix and vectors*/ + ElementMatrix* Ke = element->NewElementMatrix(MLHOApproximationEnum); + IssmDouble* basis = xNew(numnodes); + + /*Retrieve all inputs and parameters*/ + element->GetVerticesCoordinates(&xyz_list); + element->FindParam(&friction_style,GroundinglineFrictionInterpolationEnum); + Input* surface_input = element->GetInput(SurfaceEnum); _assert_(surface_input); + Input* gllevelset_input = NULL; + + /*build friction object, used later on: */ + Friction* friction=new Friction(element,dim); + + /*Recover portion of element that is grounded*/ + if(!(friction_style==SubelementFriction2Enum)) phi=element->GetGroundedPortion(xyz_list); + if(friction_style==SubelementFriction2Enum){ + gllevelset_input=element->GetInput(MaskOceanLevelsetEnum); _assert_(gllevelset_input); + element->GetGroundedPart(&point1,&fraction1,&fraction2,&mainlyfloating); + gauss = element->NewGauss(point1,fraction1,fraction2,mainlyfloating,2); + } + else{ + gauss = element->NewGauss(2); + } + /* Start looping on the number of gaussian points: */ + while(gauss->next()){ + + friction->GetAlpha2(&alpha2,gauss); + if(friction_style==SubelementFriction1Enum) alpha2=phi*alpha2; + else if(friction_style==SubelementFriction2Enum){ + gllevelset_input->GetInputValue(&gllevelset, gauss); + if(gllevelset<0.) alpha2=0.; + } + else if(friction_style==NoFrictionOnPartiallyFloatingEnum){ + if (phi<0.99999999) alpha2=0.; + } + else _error_("friction interpolation "<NodalFunctions(basis,gauss); + element->JacobianDeterminant(&Jdet, xyz_list,gauss); + + for(int i=0;ivalues[(4*i+0)*2*2*numnodes+4*j+0] += alpha2*gauss->weight*Jdet*basis[i]*basis[j]; + Ke->values[(4*i+2)*2*2*numnodes+4*j+2] += alpha2*gauss->weight*Jdet*basis[i]*basis[j]; + } + } + } + + /*Clean up and return*/ + delete gauss; + delete friction; + xDelete(xyz_list); + xDelete(basis); + return Ke; + + + //itapopo OLD - testing above /*Intermediaries*/ - int domaintype; + //int domaintype; Element* basalelement; @@ -2764,9 +2846,9 @@ //Element* basalelement = element->SpawnBasalElement(); - ElementMatrix* Ke = basalelement->NewElementMatrix(MLHOApproximationEnum); + //ElementMatrix* Ke = basalelement->NewElementMatrix(MLHOApproximationEnum); ElementMatrix* KeSSA = CreateKMatrixSSAFriction(basalelement); //only to get K11 and K33 /*Fetch number of nodes and dof for this finite element*/ - int numnodes = basalelement->GetNumberOfNodes(); + //int numnodes = basalelement->GetNumberOfNodes(); for(int i=0;iEcho(); + Ke->Echo(); + _error_("mesh "); + /*Transform Coordinate System*/ //basalelement->TransformStiffnessMatrixCoord(Ke,XYMLHOEnum); @@ -2831,12 +2917,13 @@ gauss->SynchronizeGaussBase(gauss_base); - element->JacobianDeterminant(&Jdet,xyz_list,gauss); + element->JacobianDeterminant(&Jdet,xyz_list,gauss_base); basalelement->NodalFunctionsDerivatives(dbasis,xyz_list,gauss_base); - element->NodalFunctions(basis, gauss); + basalelement->NodalFunctions(basis, gauss_base); thickness_input->GetInputValue(&thickness, gauss); n_input->GetInputValue(&n,gauss); + //FIXME testing with L1L2-type viscosity element->material->ViscosityL1L2(&viscosity,xyz_list,gauss,vx_input,vy_input,surface_input); - //viscosity=10e10;//itapopo + //viscosity=10e13;//itapopo for(int i=0;i(vx[i])) _error_("NaN found in solution vector"); if(xIsInf(vx[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"); - //if(dim==3){ - vy[i] =values[i*4+2]; //basal vy - vshy[i]=values[i*4+3]; - if(xIsNan(vy[i])) _error_("NaN found in solution vector"); - if(xIsInf(vy[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"); - //} - } - - /*Add vx and vy as inputs to the tria element (basal velocities): */ - element->AddBasalInput(VxBaseEnum,vx,element->GetElementType()); - element->AddBasalInput(VyBaseEnum,vy,element->GetElementType()); - - /*Compute suface velocities vx and vy*/ - if(domaintype==Domain2DhorizontalEnum) { - for(i=0;i(vy[i])) _error_("NaN found in solution vector"); + if(xIsInf(vy[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"); + } + + /*Add vx and vy as inputs to the tria element (shear velocities): */ + element->AddBasalInput(VxShearEnum,vshx,element->GetElementType()); + element->AddBasalInput(VyShearEnum,vshy,element->GetElementType()); + + /*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()); @@ -3191,5 +3270,5 @@ void StressbalanceAnalysis::GetSolutionFromInputsMLHO(Vector* solution,Element* element){/*{{{*/ - IssmDouble vx,vy,vbx,vby; + IssmDouble vx,vy,vshx,vshy; int domaintype,dim,approximation,dofpernode; int* doflist = NULL; @@ -3215,11 +3294,7 @@ /*Get inputs*/ Input* vx_input =element->GetInput(VxEnum); _assert_(vx_input); - Input* vxbase_input =element->GetInput(VxBaseEnum); _assert_(vxbase_input); - Input* vy_input =NULL; - Input* vybase_input =NULL; - if(domaintype!=Domain2DverticalEnum){ - vy_input =element->GetInput(VyEnum); _assert_(vy_input); - vybase_input=element->GetInput(VyBaseEnum); _assert_(vybase_input); - } + Input* vxshear_input =element->GetInput(VxShearEnum); _assert_(vxshear_input); + Input* vy_input =element->GetInput(VyEnum); _assert_(vy_input); + Input* vyshear_input =element->GetInput(VyShearEnum); _assert_(vyshear_input); /*Ok, we have vx and vy in values, fill in vx and vy arrays: */ @@ -3229,14 +3304,12 @@ /*Recover vx and vy*/ - vx_input->GetInputValue(&vx,gauss); - vxbase_input->GetInputValue(&vbx,gauss); - values[i*4+0]=vbx; //basal vx - values[i*4+1]=vx-vbx; //shear vx - //if(dofpernode==2){ - vy_input->GetInputValue(&vy,gauss); - vybase_input->GetInputValue(&vby,gauss); - values[i*4+2]=vby; //basal vy - values[i*4+3]=vy-vby; //shear vy - //} + vx_input->GetInputValue(&vx,gauss); //base vx + vxshear_input->GetInputValue(&vshx,gauss);//shear vx + values[i*4+0]=vx; //base vx + values[i*4+1]=vshx; //shear vx + vy_input->GetInputValue(&vy,gauss); //base vy + vyshear_input->GetInputValue(&vshy,gauss);//shear vy + values[i*4+2]=vy; //base vy + values[i*4+3]=vshy;//shear vy } Index: /issm/trunk-jpl/src/c/shared/Enum/Enum.vim =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/Enum.vim (revision 26191) +++ /issm/trunk-jpl/src/c/shared/Enum/Enum.vim (revision 26192) @@ -935,13 +935,13 @@ 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 VyAverageEnum -syn keyword cConstant VyBaseEnum syn keyword cConstant VyEnum syn keyword cConstant VyMeshEnum syn keyword cConstant VyObsEnum +syn keyword cConstant VyShearEnum syn keyword cConstant VzEnum syn keyword cConstant VzFSEnum @@ -1451,4 +1451,5 @@ syn keyword cType Cfsurfacesquare syn keyword cType Channel +syn keyword cType classes syn keyword cType Constraint syn keyword cType Constraints @@ -1457,6 +1458,6 @@ syn keyword cType ControlInput syn keyword cType Covertree +syn keyword cType DatasetInput syn keyword cType DataSetParam -syn keyword cType DatasetInput syn keyword cType Definition syn keyword cType DependentObject @@ -1471,6 +1472,6 @@ syn keyword cType ElementInput syn keyword cType ElementMatrix +syn keyword cType Elements syn keyword cType ElementVector -syn keyword cType Elements syn keyword cType ExponentialVariogram syn keyword cType ExternalResult @@ -1479,9 +1480,10 @@ syn keyword cType Friction syn keyword cType Gauss +syn keyword cType GaussianVariogram +syn keyword cType gaussobjects syn keyword cType GaussPenta syn keyword cType GaussSeg syn keyword cType GaussTetra syn keyword cType GaussTria -syn keyword cType GaussianVariogram syn keyword cType GenericExternalResult syn keyword cType GenericOption @@ -1498,4 +1500,5 @@ syn keyword cType IssmDirectApplicInterface syn keyword cType IssmParallelDirectApplicInterface +syn keyword cType krigingobjects syn keyword cType Load syn keyword cType Loads @@ -1508,4 +1511,5 @@ syn keyword cType Matice syn keyword cType Matlitho +syn keyword cType matrixobjects syn keyword cType MatrixParam syn keyword cType Misfit @@ -1520,6 +1524,6 @@ syn keyword cType Observations syn keyword cType Option +syn keyword cType Options syn keyword cType OptionUtilities -syn keyword cType Options syn keyword cType Param syn keyword cType Parameters @@ -1535,11 +1539,11 @@ syn keyword cType Regionaloutput syn keyword cType Results +syn keyword cType Riftfront syn keyword cType RiftStruct -syn keyword cType Riftfront syn keyword cType SealevelMasks syn keyword cType Seg syn keyword cType SegInput +syn keyword cType Segment syn keyword cType SegRef -syn keyword cType Segment syn keyword cType SpcDynamic syn keyword cType SpcStatic @@ -1560,8 +1564,4 @@ syn keyword cType Vertex syn keyword cType Vertices -syn keyword cType classes -syn keyword cType gaussobjects -syn keyword cType krigingobjects -syn keyword cType matrixobjects syn keyword cType AdjointBalancethickness2Analysis syn keyword cType AdjointBalancethicknessAnalysis Index: /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h (revision 26191) +++ /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h (revision 26192) @@ -932,13 +932,13 @@ VelEnum, VxAverageEnum, - VxBaseEnum, VxEnum, VxMeshEnum, VxObsEnum, + VxShearEnum, VyAverageEnum, - VyBaseEnum, VyEnum, VyMeshEnum, VyObsEnum, + VyShearEnum, VzEnum, VzFSEnum, Index: /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp (revision 26191) +++ /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp (revision 26192) @@ -937,13 +937,13 @@ 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 VyAverageEnum : return "VyAverage"; - case VyBaseEnum : return "VyBase"; case VyEnum : return "Vy"; case VyMeshEnum : return "VyMesh"; case VyObsEnum : return "VyObs"; + case VyShearEnum : return "VyShear"; 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 26191) +++ /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp (revision 26192) @@ -958,13 +958,13 @@ 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,"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,"Vz")==0) return VzEnum; else if (strcmp(name,"VzFS")==0) return VzFSEnum; Index: /issm/trunk-jpl/src/m/classes/flowequation.m =================================================================== --- /issm/trunk-jpl/src/m/classes/flowequation.m (revision 26191) +++ /issm/trunk-jpl/src/m/classes/flowequation.m (revision 26192) @@ -122,6 +122,6 @@ md = checkfield(md,'fieldname','flowequation.borderFS','size',[md.mesh.numberofvertices 1],'values',[0 1]); if strcmp(domaintype(md.mesh),'2Dhorizontal') - md = checkfield(md,'fieldname','flowequation.vertex_equation','size',[md.mesh.numberofvertices 1],'values',[1:2]); - md = checkfield(md,'fieldname','flowequation.element_equation','size',[md.mesh.numberofelements 1],'values',[1:2]); + md = checkfield(md,'fieldname','flowequation.vertex_equation','size',[md.mesh.numberofvertices 1],'values',[1,2,4]); + md = checkfield(md,'fieldname','flowequation.element_equation','size',[md.mesh.numberofelements 1],'values',[1,2,4]); elseif strcmp(domaintype(md.mesh),'3Dsurface') md = checkfield(md,'fieldname','flowequation.vertex_equation','size',[md.mesh.numberofvertices 1],'values',[1:2]); Index: /issm/trunk-jpl/src/m/solve/parseresultsfromdisk.m =================================================================== --- /issm/trunk-jpl/src/m/solve/parseresultsfromdisk.m (revision 26191) +++ /issm/trunk-jpl/src/m/solve/parseresultsfromdisk.m (revision 26192) @@ -224,4 +224,8 @@ field = field*yts; elseif strcmp(fieldname,'Vel'), + field = field*yts; + elseif strcmp(fieldname,'VxShear'), + field = field*yts; + elseif strcmp(fieldname,'VyShear'), field = field*yts; elseif strcmp(fieldname,'BasalforcingsGroundediceMeltingRate'),