Changeset 21620


Ignore:
Timestamp:
03/22/17 10:04:00 (8 years ago)
Author:
adhikari
Message:

CHG: compute strain rate and vorticity tensor in 2D ESA model

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

Legend:

Unmodified
Added
Removed

  • issm/trunk-jpl/src/c/classes/Elements/Element.cpp

    r21578 r21620  
    15611561                                name==EsaNmotionEnum ||
    15621562                                name==EsaEmotionEnum ||
     1563                                name==EsaStrainratexxEnum ||
     1564                                name==EsaStrainratexyEnum ||
     1565                                name==EsaStrainrateyyEnum ||
     1566                                name==EsaRotationrateEnum ||
    15631567                                name==EsaDeltathicknessEnum ||
    15641568                                name==GiaWEnum ||
    … …  
    21212125                case DeviatoricStresszzEnum:
    21222126                case DeviatoricStresseffectiveEnum: this->ComputeDeviatoricStressTensor(); break;
     2127                case EsaStrainratexxEnum:
     2128                case EsaStrainratexyEnum:
     2129                case EsaStrainrateyyEnum:
     2130                case EsaRotationrateEnum: this->ComputeEsaStrainAndVorticity(); break;
    21232131                case SigmaNNEnum: this->ComputeSigmaNN(); break;
    21242132                case LambdaSEnum: this->ComputeLambdaS(); break;
    … …  
    26732681}
    26742682/*}}}*/
     2683void       Element::StrainRateESA(IssmDouble* epsilon,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input){/*{{{*/
     2684
     2685        /*Intermediaries*/
     2686        IssmDouble dvx[3];
     2687        IssmDouble dvy[3];
     2688
     2689        /*Check that both inputs have been found*/
     2690        if(!vx_input || !vy_input){
     2691                _error_("Input missing. Here are the input pointers we have for vx: " << vx_input << ", vy: " << vy_input << "\n");
     2692        }
     2693
     2694        /*Get strain rate assuming that epsilon has been allocated*/
     2695        vx_input->GetInputDerivativeValue(&dvx[0],xyz_list,gauss);
     2696        vy_input->GetInputDerivativeValue(&dvy[0],xyz_list,gauss);
     2697        epsilon[0] = dvx[0];    // normal strain rate x-direction 
     2698        epsilon[1] = dvy[1]; // normal strain rate y-direction
     2699        epsilon[2] = 0.5*(dvx[1] + dvy[0]); // shear strain rate
     2700        epsilon[3] = 0.5*(dvx[1] - dvy[0]); // rotation rate
     2701
     2702}/*}}}*/
    26752703void       Element::StrainRateFS(IssmDouble* epsilon,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input,Input* vz_input){/*{{{*/
    26762704        /*Compute the 3d Strain Rate (6 components):

  • issm/trunk-jpl/src/c/classes/Elements/Element.h

    r21583 r21620  
    141141                int                Sid();
    142142                void               SmbGemb();
     143                void               StrainRateESA(IssmDouble* epsilon,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input);
    143144                void               StrainRateFS(IssmDouble* epsilon,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input,Input* vz_input);
    144145                void               StrainRateHO(IssmDouble* epsilon,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input);
    … …  
    182183                virtual void       ComputeSigmaNN(void)=0;
    183184                virtual void       ComputeStressTensor(void)=0;
     185                virtual void       ComputeEsaStrainAndVorticity(void)=0;
    184186                virtual void       Configure(Elements* elements,Loads* loads,Nodes* nodes,Vertices* vertices,Materials* materials,Parameters* parameters)=0;
    185187                virtual void       ControlInputSetGradient(IssmDouble* gradient,int enum_type,int control_index)=0;

  • issm/trunk-jpl/src/c/classes/Elements/Penta.h

    r21583 r21620  
    5353                void           ComputeBasalStress(Vector<IssmDouble>* sigma_b);
    5454                void           ComputeDeviatoricStressTensor();
     55                void           ComputeEsaStrainAndVorticity(){_error_("not implemented yet!");};
    5556                void           ComputeSigmaNN(){_error_("not implemented yet");};
    5657                void           ComputeStressTensor();

  • issm/trunk-jpl/src/c/classes/Elements/Seg.h

    r21583 r21620  
    4848                void        ComputeBasalStress(Vector<IssmDouble>* sigma_b){_error_("not implemented yet");};
    4949                void        ComputeDeviatoricStressTensor(){_error_("not implemented yet");};
     50                void        ComputeEsaStrainAndVorticity(){_error_("not implemented yet!");};
    5051                void        ComputeSigmaNN(){_error_("not implemented yet");};
    5152                void        ComputeStressTensor(){_error_("not implemented yet");};

  • issm/trunk-jpl/src/c/classes/Elements/Tetra.h

    r21583 r21620  
    5050                void        ComputeStressTensor(){_error_("not implemented yet");};
    5151                void        ComputeDeviatoricStressTensor(){_error_("not implemented yet");};
     52                void        ComputeEsaStrainAndVorticity(){_error_("not implemented yet!");};
    5253                void        Configure(Elements* elements,Loads* loads,Nodes* nodesin,Vertices* verticesin,Materials* materials,Parameters* parameters);
    5354                void        ControlInputSetGradient(IssmDouble* gradient,int enum_type,int control_index){_error_("not implemented yet");};

  • issm/trunk-jpl/src/c/classes/Elements/Tria.cpp

    r21583 r21620  
    410410        this->inputs->AddInput(new TriaInput(DeviatoricStresszzEnum,&tau_zz[0],P1Enum));
    411411        this->inputs->AddInput(new TriaInput(DeviatoricStresseffectiveEnum,&tau_e[0],P1Enum));
     412
     413        /*Clean up and return*/
     414        delete gauss;
     415}
     416/*}}}*/
     417void                    Tria::ComputeEsaStrainAndVorticity(){ /*{{{*/
     418
     419        IssmDouble  xyz_list[NUMVERTICES][3];
     420        IssmDouble  epsilon[4]; /* epsilon=[exx,eyy,exy+ (shear),exy- (rotation)];*/
     421        IssmDouble  strain_xx[NUMVERTICES];
     422        IssmDouble  strain_yy[NUMVERTICES];
     423        IssmDouble  strain_xy[NUMVERTICES];
     424        IssmDouble  vorticity_xy[NUMVERTICES];
     425        GaussTria*  gauss=NULL;
     426       
     427        /* Get node coordinates and dof list: */
     428        ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
     429
     430        /*Retrieve all inputs we will be needing: */
     431        Input* vx_input=this->GetInput(EsaEmotionEnum); _assert_(vx_input);
     432        Input* vy_input=this->GetInput(EsaNmotionEnum); _assert_(vy_input);
     433       
     434        /* Start looping on the number of vertices: */
     435        gauss=new GaussTria();
     436        for (int iv=0;iv<NUMVERTICES;iv++){
     437                gauss->GaussVertex(iv);
     438
     439                /*Compute strain rate and vorticity rate: */
     440                this->StrainRateESA(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input);
     441
     442                /*Compute Stress*/
     443                strain_xx[iv]=epsilon[0];
     444                strain_yy[iv]=epsilon[1];
     445                strain_xy[iv]=epsilon[2];
     446                vorticity_xy[iv]=epsilon[3];
     447        }
     448
     449        /*Add Stress tensor components into inputs*/
     450        this->inputs->AddInput(new TriaInput(EsaStrainratexxEnum,&strain_xx[0],P1Enum));
     451        this->inputs->AddInput(new TriaInput(EsaStrainrateyyEnum,&strain_yy[0],P1Enum));
     452        this->inputs->AddInput(new TriaInput(EsaStrainratexyEnum,&strain_xy[0],P1Enum));
     453        this->inputs->AddInput(new TriaInput(EsaRotationrateEnum,&vorticity_xy[0],P1Enum));
    412454
    413455        /*Clean up and return*/

  • issm/trunk-jpl/src/c/classes/Elements/Tria.h

    r21583 r21620  
    5757                void        ComputeBasalStress(Vector<IssmDouble>* sigma_b);
    5858                void        ComputeDeviatoricStressTensor();
     59                void        ComputeEsaStrainAndVorticity();
    5960                void        ComputeSigmaNN();
    6061                void        ComputeStressTensor();

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

    r21583 r21620  
    805805        EsaNmotionEnum,
    806806        EsaEmotionEnum,
     807        EsaStrainratexxEnum,
     808   EsaStrainratexyEnum,
     809   EsaStrainrateyyEnum,
     810   EsaRotationrateEnum,
    807811        EsaDeltathicknessEnum,
    808812        EsaUElasticEnum,

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

    r21583 r21620  
    783783                case EsaNmotionEnum : return "EsaNmotion";
    784784                case EsaEmotionEnum : return "EsaEmotion";
     785                case EsaStrainratexxEnum : return "EsaStrainratexx";
     786                case EsaStrainratexyEnum : return "EsaStrainratexy";
     787                case EsaStrainrateyyEnum : return "EsaStrainrateyy";
     788                case EsaRotationrateEnum : return "EsaRotationrate";
    785789                case EsaDeltathicknessEnum : return "EsaDeltathickness";
    786790                case EsaUElasticEnum : return "EsaUElastic";

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

    r21603 r21620  
    801801              else if (strcmp(name,"EsaNmotion")==0) return EsaNmotionEnum;
    802802              else if (strcmp(name,"EsaEmotion")==0) return EsaEmotionEnum;
     803              else if (strcmp(name,"EsaStrainratexx")==0) return EsaStrainratexxEnum;
     804              else if (strcmp(name,"EsaStrainratexy")==0) return EsaStrainratexyEnum;
     805              else if (strcmp(name,"EsaStrainrateyy")==0) return EsaStrainrateyyEnum;
     806              else if (strcmp(name,"EsaRotationrate")==0) return EsaRotationrateEnum;
    803807              else if (strcmp(name,"EsaDeltathickness")==0) return EsaDeltathicknessEnum;
    804808              else if (strcmp(name,"EsaUElastic")==0) return EsaUElasticEnum;
    … …  
    871875              else if (strcmp(name,"VertexPId")==0) return VertexPIdEnum;
    872876              else if (strcmp(name,"VertexSId")==0) return VertexSIdEnum;
    873               else if (strcmp(name,"Option")==0) return OptionEnum;
     877         else stage=8;
     878   }
     879   if(stage==8){
     880              if (strcmp(name,"Option")==0) return OptionEnum;
    874881              else if (strcmp(name,"GenericOption")==0) return GenericOptionEnum;
    875882              else if (strcmp(name,"OptionCell")==0) return OptionCellEnum;
    876883              else if (strcmp(name,"OptionStruct")==0) return OptionStructEnum;
    877          else stage=8;
    878    }
    879    if(stage==8){
    880               if (strcmp(name,"AdjointBalancethicknessAnalysis")==0) return AdjointBalancethicknessAnalysisEnum;
     884              else if (strcmp(name,"AdjointBalancethicknessAnalysis")==0) return AdjointBalancethicknessAnalysisEnum;
    881885              else if (strcmp(name,"AdjointBalancethickness2Analysis")==0) return AdjointBalancethickness2AnalysisEnum;
    882886              else if (strcmp(name,"AdjointHorizAnalysis")==0) return AdjointHorizAnalysisEnum;
    … …  
    994998              else if (strcmp(name,"Open")==0) return OpenEnum;
    995999              else if (strcmp(name,"Air")==0) return AirEnum;
    996               else if (strcmp(name,"Ice")==0) return IceEnum;
     1000         else stage=9;
     1001   }
     1002   if(stage==9){
     1003              if (strcmp(name,"Ice")==0) return IceEnum;
    9971004              else if (strcmp(name,"Melange")==0) return MelangeEnum;
    9981005              else if (strcmp(name,"Water")==0) return WaterEnum;
    9991006              else if (strcmp(name,"DataSet")==0) return DataSetEnum;
    1000          else stage=9;
    1001    }
    1002    if(stage==9){
    1003               if (strcmp(name,"Constraints")==0) return ConstraintsEnum;
     1007              else if (strcmp(name,"Constraints")==0) return ConstraintsEnum;
    10041008              else if (strcmp(name,"Loads")==0) return LoadsEnum;
    10051009              else if (strcmp(name,"Materials")==0) return MaterialsEnum;

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

    r21260 r21620  
    6868                        fielddisplay(self,'degacc','accuracy (default .01 deg) for numerical discretization of the Green''s functions');
    6969                        fielddisplay(self,'transitions','indices into parts of the mesh that will be icecaps');
    70                         fielddisplay(self,'requested_outputs','additional outputs requested (default: EsaUmotion)');
     70                        fielddisplay(self,'requested_outputs','additional outputs requested (e.g., EsaUmotion, EsaStrainratexx, EsaRotationrate)');
    7171
    7272                end % }}}
Note: See TracChangeset for help on using the changeset viewer.