Changeset 26192


Ignore:
Timestamp:
04/16/21 08:52:05 (4 years ago)
Author:
tsantos
Message:

CHG: working on Mono Layer. Keeping Shear velocities in inputs (VxShear/VyShear). VxBase/VyBase are now Vx/Vy

Location:
issm/trunk-jpl/src
Files:
7 edited

Legend:

Unmodified
Added
Removed

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

    r26187 r26192  
    194194                                }
    195195                        }
    196                         else{//itapopo testing here
     196                        else{//itapopo FIXME only base velocities are constrained here
    197197                                IoModelToConstraintsx(constraints,iomodel,"md.stressbalance.spcvx",StressbalanceAnalysisEnum,finiteelement,0);
    198                                 IoModelToConstraintsx(constraints,iomodel,"md.stressbalance.spcvx",StressbalanceAnalysisEnum,finiteelement,1);
     198                                //IoModelToConstraintsx(constraints,iomodel,"md.stressbalance.spcvx",StressbalanceAnalysisEnum,finiteelement,1);
    199199                                if(iomodel->domaintype!=Domain2DverticalEnum){
    200200                                        IoModelToConstraintsx(constraints,iomodel,"md.stressbalance.spcvy",StressbalanceAnalysisEnum,finiteelement,2);
    201                                         IoModelToConstraintsx(constraints,iomodel,"md.stressbalance.spcvy",StressbalanceAnalysisEnum,finiteelement,3);
     201                                        //IoModelToConstraintsx(constraints,iomodel,"md.stressbalance.spcvy",StressbalanceAnalysisEnum,finiteelement,3);
    202202                                }
    203203                        }
    … …  
    783783        iomodel->FetchDataToInput(inputs,elements,"md.initialization.vx",VxEnum,0.);
    784784        iomodel->FetchDataToInput(inputs,elements,"md.initialization.vy",VyEnum,0.);
    785         if(isMLHO){//itapopo
    786                 iomodel->FetchDataToInput(inputs,elements,"md.initialization.vx",VxBaseEnum,0.);
    787                 iomodel->FetchDataToInput(inputs,elements,"md.initialization.vy",VyBaseEnum,0.);
     785        if(isMLHO){//itapopo FIXME the shear velocities should be initialized correctly
     786                iomodel->FetchDataToInput(inputs,elements,"md.initialization.vx",VxShearEnum,0.);
     787                iomodel->FetchDataToInput(inputs,elements,"md.initialization.vy",VyShearEnum,0.);
    788788        }
    789789        iomodel->FetchDataToInput(inputs,elements,"md.stressbalance.loadingforcex",LoadingforceXEnum);
    … …  
    20842084        /*Use the dof list to index into the solution vector: */
    20852085        for(i=0;i<numdof;i++) values[i]=solution[doflist[i]];
    2086 
     2086       
    20872087        /*Transform solution in Cartesian Space*/
    20882088        if(dim==2) basalelement->TransformSolutionCoord(&values[0],XYEnum);
    … …  
    27292729/*MLHO*/
    27302730ElementMatrix* StressbalanceAnalysis::CreateKMatrixMLHO(Element* element){/*{{{*/
    2731 
     2731       
     2732        _error_("Mono Layer Higher-Order called, not fully tested. If you are sure in using it, comment this line.");
     2733       
    27322734        /* Check if ice in element */
    27332735        if(!element->IsIceInElement()) return NULL;
    … …  
    27452747
    27462748        if(!element->IsOnBase() || element->IsFloating()) return NULL;
     2749
     2750        /*Intermediaries*/
     2751        int         dim,domaintype;
     2752        bool        mainlyfloating;
     2753        int         friction_style,point1;
     2754        IssmDouble  alpha2,Jdet,fraction1,fraction2;
     2755        IssmDouble  gllevelset,phi=1.;
     2756        IssmDouble *xyz_list  = NULL;
     2757        Gauss*      gauss     = NULL;
     2758
     2759        /*Get problem dimension*/
     2760        element->FindParam(&domaintype,DomainTypeEnum);
     2761        switch(domaintype){
     2762                case Domain2DverticalEnum:   dim = 1;break;
     2763                case Domain2DhorizontalEnum: dim = 2;break;
     2764                case Domain3DEnum:           dim = 2;break;
     2765                default: _error_("mesh "<<EnumToStringx(domaintype)<<" not supported yet");
     2766        }
     2767
     2768        /*Fetch number of nodes and dof for this finite element*/
     2769        int numnodes = element->GetNumberOfNodes();
     2770        int numdof   = numnodes*dim;
     2771
     2772        /*Initialize Element matrix and vectors*/
     2773        ElementMatrix* Ke = element->NewElementMatrix(MLHOApproximationEnum);
     2774        IssmDouble*    basis  = xNew<IssmDouble>(numnodes);
     2775
     2776        /*Retrieve all inputs and parameters*/
     2777        element->GetVerticesCoordinates(&xyz_list);
     2778        element->FindParam(&friction_style,GroundinglineFrictionInterpolationEnum);
     2779        Input* surface_input    = element->GetInput(SurfaceEnum); _assert_(surface_input);
     2780        Input* gllevelset_input = NULL;
     2781
     2782        /*build friction object, used later on: */
     2783        Friction* friction=new Friction(element,dim);
     2784
     2785        /*Recover portion of element that is grounded*/
     2786        if(!(friction_style==SubelementFriction2Enum)) phi=element->GetGroundedPortion(xyz_list);
     2787        if(friction_style==SubelementFriction2Enum){
     2788                gllevelset_input=element->GetInput(MaskOceanLevelsetEnum); _assert_(gllevelset_input);
     2789                element->GetGroundedPart(&point1,&fraction1,&fraction2,&mainlyfloating);
     2790           gauss = element->NewGauss(point1,fraction1,fraction2,mainlyfloating,2);
     2791        }
     2792        else{
     2793                gauss = element->NewGauss(2);
     2794        }
    27472795       
     2796        /* Start  looping on the number of gaussian points: */
     2797        while(gauss->next()){
     2798
     2799                friction->GetAlpha2(&alpha2,gauss);
     2800                if(friction_style==SubelementFriction1Enum) alpha2=phi*alpha2;
     2801                else if(friction_style==SubelementFriction2Enum){
     2802                        gllevelset_input->GetInputValue(&gllevelset, gauss);
     2803                        if(gllevelset<0.) alpha2=0.;
     2804                }
     2805                else if(friction_style==NoFrictionOnPartiallyFloatingEnum){
     2806                        if (phi<0.99999999) alpha2=0.;
     2807                }
     2808                else  _error_("friction interpolation "<<EnumToStringx(friction_style)<<" not implemented yet");
     2809
     2810                element->NodalFunctions(basis,gauss);
     2811                element->JacobianDeterminant(&Jdet, xyz_list,gauss);
     2812
     2813                for(int i=0;i<numnodes;i++){
     2814                        for(int j=0;j<numnodes;j++){
     2815                                Ke->values[(4*i+0)*2*2*numnodes+4*j+0] += alpha2*gauss->weight*Jdet*basis[i]*basis[j];
     2816                                Ke->values[(4*i+2)*2*2*numnodes+4*j+2] += alpha2*gauss->weight*Jdet*basis[i]*basis[j];
     2817                        }
     2818                }
     2819        }
     2820
     2821        /*Clean up and return*/
     2822        delete gauss;
     2823        delete friction;
     2824        xDelete<IssmDouble>(xyz_list);
     2825        xDelete<IssmDouble>(basis);
     2826        return Ke;
     2827
     2828
     2829        //itapopo OLD - testing above
    27482830        /*Intermediaries*/
    2749         int      domaintype;
     2831        //int      domaintype;
    27502832   Element* basalelement;
    27512833
    … …  
    27642846
    27652847        //Element* basalelement = element->SpawnBasalElement();
    2766         ElementMatrix* Ke    = basalelement->NewElementMatrix(MLHOApproximationEnum);
     2848        //ElementMatrix* Ke    = basalelement->NewElementMatrix(MLHOApproximationEnum);
    27672849        ElementMatrix* KeSSA = CreateKMatrixSSAFriction(basalelement); //only to get K11 and K33
    27682850       
    27692851        /*Fetch number of nodes and dof for this finite element*/
    2770         int numnodes = basalelement->GetNumberOfNodes();
     2852        //int numnodes = basalelement->GetNumberOfNodes();
    27712853
    27722854        for(int i=0;i<numnodes;i++){
    … …  
    27762858      } 
    27772859   }
    2778        
     2860
     2861        KeSSA->Echo();
     2862        Ke->Echo();
     2863        _error_("mesh ");
     2864
    27792865        /*Transform Coordinate System*/
    27802866        //basalelement->TransformStiffnessMatrixCoord(Ke,XYMLHOEnum);
    … …  
    28312917                gauss->SynchronizeGaussBase(gauss_base);
    28322918
    2833                 element->JacobianDeterminant(&Jdet,xyz_list,gauss);
     2919                element->JacobianDeterminant(&Jdet,xyz_list,gauss_base);
    28342920                basalelement->NodalFunctionsDerivatives(dbasis,xyz_list,gauss_base);
    2835                 element->NodalFunctions(basis, gauss);
     2921                basalelement->NodalFunctions(basis, gauss_base);
    28362922
    28372923                thickness_input->GetInputValue(&thickness, gauss);
    28382924                n_input->GetInputValue(&n,gauss);
     2925                //FIXME testing with L1L2-type viscosity
    28392926                element->material->ViscosityL1L2(&viscosity,xyz_list,gauss,vx_input,vy_input,surface_input);
    2840                 //viscosity=10e10;//itapopo
     2927                //viscosity=10e13;//itapopo
    28412928
    28422929                for(int i=0;i<numnodes;i++){//shape functions on tria element
    … …  
    30403127                base_under_water    = min(0.,base-sealevel);           // 0 if the bottom of the glacier is above water level
    30413128                /*Vertically integrated pressure - SSA type*/
    3042                 water_pressure = 1.0/2.0*gravity*rho_water*(surface_under_water*surface_under_water - base_under_water*base_under_water);
    3043                 ice_pressure   = 1.0/2.0*gravity*rho_ice*thickness*thickness;
     3129                water_pressure = (1.0/2.0)*gravity*rho_water*(surface_under_water*surface_under_water - base_under_water*base_under_water);
     3130                ice_pressure   = (1.0/2.0)*gravity*rho_ice*thickness*thickness;
    30443131                pressure = ice_pressure + water_pressure;
    30453132                /*Vertically integrated pressure - HO type*/
    30463133                b=min(0.,base-sealevel); // 0 if the bottom of the glacier is above water level
    30473134                s=thickness+b; // ice surface regardless of whether the top of the glacier is above water level or not
    3048                 water_pressure_sh = gravity*rho_water*(-b*b/2 + pow(-b,n+2)*( -s/(n+2) -b/(n+3) )/pow(thickness,n+1));
     3135                water_pressure_sh = gravity*rho_water*( -b*b/2 + (pow(-b,n+2)/pow(thickness,n+1))*(-s/(n+2)-b/(n+3)) );
    30493136                ice_pressure_sh   = gravity*rho_ice*thickness*thickness*(n+1)/(2*(n+3));
    30503137                pressure_sh = ice_pressure_sh + water_pressure_sh;
    … …  
    31313218        for(i=0;i<numdof;i++) values[i]=solution[doflist[i]];
    31323219
    3133         std::cout<<"********  Element ID="<<element->Id()<<"\n";
    3134         for(i=0;i<numdof;i++){
    3135                 std::cout<<values[i]<<"\n";
    3136         }
     3220        //std::cout<<"MLHO  ********  Element ID="<<element->Id()<<"\n";
     3221        //for(i=0;i<numdof;i++){
     3222        //      std::cout<<values[i]*31536000<<"\n";
     3223        //}
    31373224
    31383225        /*Transform solution in Cartesian Space*/
    … …  
    31423229   /*Ok, we have vx and vy in values, fill in vx and vy arrays: */
    31433230   for(i=0;i<numnodes;i++){ //numnodes of the 2D mesh in which the MLHO is written
    3144       vx[i]  =values[i*4+0]; //basal vx
    3145       vshx[i]=values[i*4+1];
     3231      vx[i]  =values[i*4+0]; //base vx
     3232      vshx[i]=values[i*4+1]; //shear vx
    31463233                if(xIsNan<IssmDouble>(vx[i]))           _error_("NaN found in solution vector");
    31473234      if(xIsInf<IssmDouble>(vx[i]))             _error_("Inf found in solution vector");
    31483235                if(xIsNan<IssmDouble>(vshx[i])) _error_("NaN found in solution vector");
    31493236      if(xIsInf<IssmDouble>(vshx[i]))   _error_("Inf found in solution vector");
    3150       //if(dim==3){
    3151          vy[i] =values[i*4+2]; //basal vy
    3152          vshy[i]=values[i*4+3];
    3153          if(xIsNan<IssmDouble>(vy[i]))          _error_("NaN found in solution vector");
    3154          if(xIsInf<IssmDouble>(vy[i]))          _error_("Inf found in solution vector");
    3155          if(xIsNan<IssmDouble>(vshy[i]))        _error_("NaN found in solution vector");
    3156          if(xIsInf<IssmDouble>(vshy[i]))        _error_("Inf found in solution vector");
    3157       //}
    3158    }
    3159 
    3160         /*Add vx and vy as inputs to the tria element (basal velocities): */
    3161         element->AddBasalInput(VxBaseEnum,vx,element->GetElementType());
    3162         element->AddBasalInput(VyBaseEnum,vy,element->GetElementType());
    3163 
    3164         /*Compute suface velocities vx and vy*/
    3165         if(domaintype==Domain2DhorizontalEnum) {
    3166                 for(i=0;i<numnodes;i++) vx[i]=vx[i]+vshx[i];
    3167                 for(i=0;i<numnodes;i++) vy[i]=vy[i]+vshy[i];
    3168         }
    3169 
    3170         /*Get Vz and compute vel (surface)*/
     3237                vy[i]  =values[i*4+2]; //base vy
     3238                vshy[i]=values[i*4+3]; //shear vy
     3239                if(xIsNan<IssmDouble>(vy[i]))           _error_("NaN found in solution vector");
     3240                if(xIsInf<IssmDouble>(vy[i]))           _error_("Inf found in solution vector");
     3241                if(xIsNan<IssmDouble>(vshy[i])) _error_("NaN found in solution vector");
     3242                if(xIsInf<IssmDouble>(vshy[i])) _error_("Inf found in solution vector");
     3243        }
     3244
     3245        /*Add vx and vy as inputs to the tria element (shear velocities): */
     3246        element->AddBasalInput(VxShearEnum,vshx,element->GetElementType());
     3247        element->AddBasalInput(VyShearEnum,vshy,element->GetElementType());
     3248
     3249        /*Get Vz and compute vel (base)*/
    31713250        basalelement->GetInputListOnNodes(&vz[0],VzEnum,0.);
    31723251        for(i=0;i<numnodes;i++) vel[i]=sqrt(vx[i]*vx[i] + vy[i]*vy[i] + vz[i]*vz[i]);
    31733252
    3174         /*Add vx and vy as inputs to the tria element surface velocities): */
     3253        /*Add vx and vy as inputs to the tria element (base velocities): */
    31753254        element->AddBasalInput(VxEnum,vx,element->GetElementType());
    31763255        element->AddBasalInput(VyEnum,vy,element->GetElementType());
    … …  
    31913270void           StressbalanceAnalysis::GetSolutionFromInputsMLHO(Vector<IssmDouble>* solution,Element* element){/*{{{*/
    31923271
    3193         IssmDouble   vx,vy,vbx,vby;
     3272        IssmDouble   vx,vy,vshx,vshy;
    31943273        int          domaintype,dim,approximation,dofpernode;
    31953274        int*         doflist = NULL;
    … …  
    32153294        /*Get inputs*/
    32163295        Input* vx_input         =element->GetInput(VxEnum);             _assert_(vx_input);
    3217         Input* vxbase_input     =element->GetInput(VxBaseEnum); _assert_(vxbase_input);
    3218         Input* vy_input         =NULL;
    3219         Input* vybase_input     =NULL;
    3220         if(domaintype!=Domain2DverticalEnum){
    3221                 vy_input                =element->GetInput(VyEnum);             _assert_(vy_input);
    3222                 vybase_input=element->GetInput(VyBaseEnum);     _assert_(vybase_input);
    3223         }
     3296        Input* vxshear_input    =element->GetInput(VxShearEnum);        _assert_(vxshear_input);
     3297        Input* vy_input         =element->GetInput(VyEnum);             _assert_(vy_input);
     3298        Input* vyshear_input    =element->GetInput(VyShearEnum);        _assert_(vyshear_input);
    32243299
    32253300        /*Ok, we have vx and vy in values, fill in vx and vy arrays: */
    … …  
    32293304
    32303305                /*Recover vx and vy*/
    3231                 vx_input->GetInputValue(&vx,gauss);
    3232                 vxbase_input->GetInputValue(&vbx,gauss);
    3233                 values[i*4+0]=vbx; //basal vx
    3234                 values[i*4+1]=vx-vbx; //shear vx
    3235                 //if(dofpernode==2){
    3236                         vy_input->GetInputValue(&vy,gauss);
    3237                         vybase_input->GetInputValue(&vby,gauss);
    3238                         values[i*4+2]=vby; //basal vy
    3239                         values[i*4+3]=vy-vby; //shear vy 
    3240                 //}
     3306                vx_input->GetInputValue(&vx,gauss);                     //base vx
     3307                vxshear_input->GetInputValue(&vshx,gauss);//shear vx
     3308                values[i*4+0]=vx;   //base vx
     3309                values[i*4+1]=vshx; //shear vx
     3310                vy_input->GetInputValue(&vy,gauss);                     //base vy
     3311                vyshear_input->GetInputValue(&vshy,gauss);//shear vy
     3312                values[i*4+2]=vy;               //base vy
     3313                values[i*4+3]=vshy;//shear vy 
    32413314        }
    32423315

  • issm/trunk-jpl/src/c/shared/Enum/Enum.vim

    r26165 r26192  
    935935syn keyword cConstant VelEnum
    936936syn keyword cConstant VxAverageEnum
    937 syn keyword cConstant VxBaseEnum
    938937syn keyword cConstant VxEnum
    939938syn keyword cConstant VxMeshEnum
    940939syn keyword cConstant VxObsEnum
     940syn keyword cConstant VxShearEnum
    941941syn keyword cConstant VyAverageEnum
    942 syn keyword cConstant VyBaseEnum
    943942syn keyword cConstant VyEnum
    944943syn keyword cConstant VyMeshEnum
    945944syn keyword cConstant VyObsEnum
     945syn keyword cConstant VyShearEnum
    946946syn keyword cConstant VzEnum
    947947syn keyword cConstant VzFSEnum
    … …  
    14511451syn keyword cType Cfsurfacesquare
    14521452syn keyword cType Channel
     1453syn keyword cType classes
    14531454syn keyword cType Constraint
    14541455syn keyword cType Constraints
    … …  
    14571458syn keyword cType ControlInput
    14581459syn keyword cType Covertree
     1460syn keyword cType DatasetInput
    14591461syn keyword cType DataSetParam
    1460 syn keyword cType DatasetInput
    14611462syn keyword cType Definition
    14621463syn keyword cType DependentObject
    … …  
    14711472syn keyword cType ElementInput
    14721473syn keyword cType ElementMatrix
     1474syn keyword cType Elements
    14731475syn keyword cType ElementVector
    1474 syn keyword cType Elements
    14751476syn keyword cType ExponentialVariogram
    14761477syn keyword cType ExternalResult
    … …  
    14791480syn keyword cType Friction
    14801481syn keyword cType Gauss
     1482syn keyword cType GaussianVariogram
     1483syn keyword cType gaussobjects
    14811484syn keyword cType GaussPenta
    14821485syn keyword cType GaussSeg
    14831486syn keyword cType GaussTetra
    14841487syn keyword cType GaussTria
    1485 syn keyword cType GaussianVariogram
    14861488syn keyword cType GenericExternalResult
    14871489syn keyword cType GenericOption
    … …  
    14981500syn keyword cType IssmDirectApplicInterface
    14991501syn keyword cType IssmParallelDirectApplicInterface
     1502syn keyword cType krigingobjects
    15001503syn keyword cType Load
    15011504syn keyword cType Loads
    … …  
    15081511syn keyword cType Matice
    15091512syn keyword cType Matlitho
     1513syn keyword cType matrixobjects
    15101514syn keyword cType MatrixParam
    15111515syn keyword cType Misfit
    … …  
    15201524syn keyword cType Observations
    15211525syn keyword cType Option
     1526syn keyword cType Options
    15221527syn keyword cType OptionUtilities
    1523 syn keyword cType Options
    15241528syn keyword cType Param
    15251529syn keyword cType Parameters
    … …  
    15351539syn keyword cType Regionaloutput
    15361540syn keyword cType Results
     1541syn keyword cType Riftfront
    15371542syn keyword cType RiftStruct
    1538 syn keyword cType Riftfront
    15391543syn keyword cType SealevelMasks
    15401544syn keyword cType Seg
    15411545syn keyword cType SegInput
     1546syn keyword cType Segment
    15421547syn keyword cType SegRef
    1543 syn keyword cType Segment
    15441548syn keyword cType SpcDynamic
    15451549syn keyword cType SpcStatic
    … …  
    15601564syn keyword cType Vertex
    15611565syn keyword cType Vertices
    1562 syn keyword cType classes
    1563 syn keyword cType gaussobjects
    1564 syn keyword cType krigingobjects
    1565 syn keyword cType matrixobjects
    15661566syn keyword cType AdjointBalancethickness2Analysis
    15671567syn keyword cType AdjointBalancethicknessAnalysis

  • issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h

    r26165 r26192  
    932932        VelEnum,
    933933        VxAverageEnum,
    934         VxBaseEnum,
    935934        VxEnum,
    936935        VxMeshEnum,
    937936        VxObsEnum,
     937        VxShearEnum,
    938938        VyAverageEnum,
    939         VyBaseEnum,
    940939        VyEnum,
    941940        VyMeshEnum,
    942941        VyObsEnum,
     942        VyShearEnum,
    943943        VzEnum,
    944944        VzFSEnum,

  • issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp

    r26165 r26192  
    937937                case VelEnum : return "Vel";
    938938                case VxAverageEnum : return "VxAverage";
    939                 case VxBaseEnum : return "VxBase";
    940939                case VxEnum : return "Vx";
    941940                case VxMeshEnum : return "VxMesh";
    942941                case VxObsEnum : return "VxObs";
     942                case VxShearEnum : return "VxShear";
    943943                case VyAverageEnum : return "VyAverage";
    944                 case VyBaseEnum : return "VyBase";
    945944                case VyEnum : return "Vy";
    946945                case VyMeshEnum : return "VyMesh";
    947946                case VyObsEnum : return "VyObs";
     947                case VyShearEnum : return "VyShear";
    948948                case VzEnum : return "Vz";
    949949                case VzFSEnum : return "VzFS";

  • issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp

    r26165 r26192  
    958958              else if (strcmp(name,"Vel")==0) return VelEnum;
    959959              else if (strcmp(name,"VxAverage")==0) return VxAverageEnum;
    960               else if (strcmp(name,"VxBase")==0) return VxBaseEnum;
    961960              else if (strcmp(name,"Vx")==0) return VxEnum;
    962961              else if (strcmp(name,"VxMesh")==0) return VxMeshEnum;
    963962              else if (strcmp(name,"VxObs")==0) return VxObsEnum;
     963              else if (strcmp(name,"VxShear")==0) return VxShearEnum;
    964964              else if (strcmp(name,"VyAverage")==0) return VyAverageEnum;
    965               else if (strcmp(name,"VyBase")==0) return VyBaseEnum;
    966965              else if (strcmp(name,"Vy")==0) return VyEnum;
    967966              else if (strcmp(name,"VyMesh")==0) return VyMeshEnum;
    968967              else if (strcmp(name,"VyObs")==0) return VyObsEnum;
     968              else if (strcmp(name,"VyShear")==0) return VyShearEnum;
    969969              else if (strcmp(name,"Vz")==0) return VzEnum;
    970970              else if (strcmp(name,"VzFS")==0) return VzFSEnum;

  • issm/trunk-jpl/src/m/classes/flowequation.m

    r25688 r26192  
    122122                        md = checkfield(md,'fieldname','flowequation.borderFS','size',[md.mesh.numberofvertices 1],'values',[0 1]);
    123123                        if strcmp(domaintype(md.mesh),'2Dhorizontal')
    124                                 md = checkfield(md,'fieldname','flowequation.vertex_equation','size',[md.mesh.numberofvertices 1],'values',[1:2]);
    125                                 md = checkfield(md,'fieldname','flowequation.element_equation','size',[md.mesh.numberofelements 1],'values',[1:2]);
     124                                md = checkfield(md,'fieldname','flowequation.vertex_equation','size',[md.mesh.numberofvertices 1],'values',[1,2,4]);
     125                                md = checkfield(md,'fieldname','flowequation.element_equation','size',[md.mesh.numberofelements 1],'values',[1,2,4]);
    126126                        elseif strcmp(domaintype(md.mesh),'3Dsurface')
    127127                                md = checkfield(md,'fieldname','flowequation.vertex_equation','size',[md.mesh.numberofvertices 1],'values',[1:2]);

  • issm/trunk-jpl/src/m/solve/parseresultsfromdisk.m

    r25714 r26192  
    224224                field = field*yts;
    225225        elseif strcmp(fieldname,'Vel'),
     226                field = field*yts;
     227        elseif strcmp(fieldname,'VxShear'),
     228                field = field*yts;
     229        elseif strcmp(fieldname,'VyShear'),
    226230                field = field*yts;
    227231        elseif strcmp(fieldname,'BasalforcingsGroundediceMeltingRate'),
Note: See TracChangeset for help on using the changeset viewer.