Changeset 17680


Ignore:
Timestamp:
04/08/14 17:20:10 (11 years ago)
Author:
cborstad
Message:

CHG: do not compute stress tensor for a 2d horiz mesh, new function ction to compute deviatoric stress

Location:
issm/trunk-jpl/src/c
Files:
8 edited

Legend:

Unmodified
Added
Removed

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

    r17674 r17680  
    430430        IssmDouble s_xx,s_xy,s_yy,s1,s2,stmp;
    431431        IssmDouble J2s,Chi,Psi,PosPsi,NegPsi;
    432         IssmDouble damage,sigma_xx,sigma_xy,sigma_yy;
     432        IssmDouble damage,tau_xx,tau_xy,tau_yy;
    433433        int equivstress;
    434434
    … …  
    445445
    446446        /*Compute stress tensor: */
    447         element->ComputeStressTensor();
     447        element->ComputeDeviatoricStressTensor();
    448448
    449449        /*retrieve what we need: */
    450         Input* sigma_xx_input  = element->GetInput(StressTensorxxEnum);     _assert_(sigma_xx_input);
    451         Input* sigma_xy_input  = element->GetInput(StressTensorxyEnum);     _assert_(sigma_xy_input);
    452         Input* sigma_yy_input  = element->GetInput(StressTensoryyEnum);     _assert_(sigma_yy_input);
     450        Input* tau_xx_input  = element->GetInput(DeviatoricStressxxEnum);     _assert_(tau_xx_input);
     451        Input* tau_xy_input  = element->GetInput(DeviatoricStressxyEnum);     _assert_(tau_xy_input);
     452        Input* tau_yy_input  = element->GetInput(DeviatoricStressyyEnum);     _assert_(tau_yy_input);
    453453        Input* damage_input    = element->GetInput(DamageDEnum); _assert_(damage_input);
    454454
    … …  
    462462               
    463463                damage_input->GetInputValue(&damage,gauss);
    464                 sigma_xx_input->GetInputValue(&sigma_xx,gauss);
    465                 sigma_xy_input->GetInputValue(&sigma_xy,gauss);
    466                 sigma_yy_input->GetInputValue(&sigma_yy,gauss);
     464                tau_xx_input->GetInputValue(&tau_xx,gauss);
     465                tau_xy_input->GetInputValue(&tau_xy,gauss);
     466                tau_yy_input->GetInputValue(&tau_yy,gauss);
    467467       
    468468                /*Calculate effective stress components*/
    469                 s_xx=sigma_xx/(1.-damage);
    470                 s_xy=sigma_xy/(1.-damage);
    471                 s_yy=sigma_yy/(1.-damage);
     469                s_xx=tau_xx/(1.-damage);
     470                s_xy=tau_xy/(1.-damage);
     471                s_yy=tau_yy/(1.-damage);
    472472
    473473                /*Calculate principal effective stresses*/

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

    r17671 r17680  
    203203                virtual void   ComputeStrainRate(Vector<IssmDouble>* eps)=0;
    204204                virtual void   ComputeStressTensor(void)=0;
     205                virtual void   ComputeDeviatoricStressTensor(void)=0;
    205206
    206207                virtual void   Update(int index, IoModel* iomodel,int analysis_counter,int analysis_type,int finite_element)=0;

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

    r17671 r17680  
    332332        this->inputs->AddInput(new PentaInput(StressTensoryzEnum,&sigma_yz[0],P1Enum));
    333333        this->inputs->AddInput(new PentaInput(StressTensorzzEnum,&sigma_zz[0],P1Enum));
     334
     335        /*Clean up and return*/
     336        delete gauss;
     337}
     338/*}}}*/
     339/*FUNCTION Penta::ComputeDeviatoricStressTensor {{{*/
     340void  Penta::ComputeDeviatoricStressTensor(){
     341
     342        IssmDouble      xyz_list[NUMVERTICES][3];
     343        IssmDouble      viscosity;
     344        IssmDouble      epsilon[6]; /* epsilon=[exx,eyy,exy];*/
     345        IssmDouble      tau_xx[NUMVERTICES];
     346        IssmDouble              tau_yy[NUMVERTICES];
     347        IssmDouble              tau_zz[NUMVERTICES];
     348        IssmDouble      tau_xy[NUMVERTICES];
     349        IssmDouble              tau_xz[NUMVERTICES];
     350        IssmDouble              tau_yz[NUMVERTICES];
     351        GaussPenta* gauss=NULL;
     352
     353        /* Get node coordinates and dof list: */
     354        ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
     355
     356        /*Retrieve all inputs we will be needing: */
     357        Input* vx_input=inputs->GetInput(VxEnum);             _assert_(vx_input);
     358        Input* vy_input=inputs->GetInput(VyEnum);             _assert_(vy_input);
     359        Input* vz_input=inputs->GetInput(VzEnum);             _assert_(vz_input);
     360
     361        /* Start looping on the number of vertices: */
     362        gauss=new GaussPenta();
     363        for (int iv=0;iv<NUMVERTICES;iv++){
     364                gauss->GaussVertex(iv);
     365
     366                /*Compute strain rate viscosity and pressure: */
     367                this->StrainRateFS(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input,vz_input);
     368                this->ViscosityFS(&viscosity,3,&xyz_list[0][0],gauss,vx_input,vy_input,vz_input);
     369
     370                /*Compute Stress*/
     371                tau_xx[iv]=2*viscosity*epsilon[0]; // tau = nu eps
     372                tau_yy[iv]=2*viscosity*epsilon[1];
     373                tau_zz[iv]=2*viscosity*epsilon[2];
     374                tau_xy[iv]=2*viscosity*epsilon[3];
     375                tau_xz[iv]=2*viscosity*epsilon[4];
     376                tau_yz[iv]=2*viscosity*epsilon[5];
     377        }
     378
     379        /*Add Stress tensor components into inputs*/
     380        this->inputs->AddInput(new PentaInput(StressTensorxxEnum,&tau_xx[0],P1Enum));
     381        this->inputs->AddInput(new PentaInput(StressTensorxyEnum,&tau_xy[0],P1Enum));
     382        this->inputs->AddInput(new PentaInput(StressTensorxzEnum,&tau_xz[0],P1Enum));
     383        this->inputs->AddInput(new PentaInput(StressTensoryyEnum,&tau_yy[0],P1Enum));
     384        this->inputs->AddInput(new PentaInput(StressTensoryzEnum,&tau_yz[0],P1Enum));
     385        this->inputs->AddInput(new PentaInput(StressTensorzzEnum,&tau_zz[0],P1Enum));
    334386
    335387        /*Clean up and return*/

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

    r17671 r17680  
    5757                void   ComputeStrainRate(Vector<IssmDouble>* eps);
    5858                void   ComputeStressTensor();
     59                void   ComputeDeviatoricStressTensor();
    5960                void   Configure(Elements* elements,Loads* loads,Nodes* nodes,Vertices* vertices,Materials* materials,Parameters* parameters);
    6061                void   ElementSizes(IssmDouble* hx,IssmDouble* hy,IssmDouble* hz);

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

    r17671 r17680  
    5656                void        ComputeStrainRate(Vector<IssmDouble>* eps){_error_("not implemented yet");};
    5757                void        ComputeStressTensor(){_error_("not implemented yet");};
     58                void        ComputeDeviatoricStressTensor(){_error_("not implemented yet");};
    5859                void        Configure(Elements* elements,Loads* loads,Nodes* nodesin,Vertices* verticesin,Materials* materials,Parameters* parameters){_error_("not implemented yet");};
    5960                void        SetCurrentConfiguration(Elements* elements,Loads* loads,Nodes* nodes,Materials* materials,Parameters* parameters){_error_("not implemented yet");};

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

    r17671 r17680  
    5656                void        ComputeStrainRate(Vector<IssmDouble>* eps){_error_("not implemented yet");};
    5757                void        ComputeStressTensor(){_error_("not implemented yet");};
     58                void        ComputeDeviatoricStressTensor(){_error_("not implemented yet");};
    5859                void        Configure(Elements* elements,Loads* loads,Nodes* nodesin,Vertices* verticesin,Materials* materials,Parameters* parameters);
    5960                void        SetCurrentConfiguration(Elements* elements,Loads* loads,Nodes* nodes,Materials* materials,Parameters* parameters);

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

    r17674 r17680  
    165165        IssmDouble      sigma_yz[NUMVERTICES]={0,0,0};
    166166        GaussTria*  gauss=NULL;
    167         int dim=2;
     167        int meshtype,dim=2;
    168168
    169169        /* Get node coordinates and dof list: */
    … …  
    171171
    172172        /*Retrieve all inputs we will be needing: */
     173        this->FindParam(&meshtype,MeshTypeEnum);
     174        if(meshtype==Mesh2DhorizontalEnum) _error_("stress tensor calculation not supported for mesh of type " <<EnumToStringx(meshtype)<<", extrude mesh or call ComputeDeviatoricStressTensor");
    173175        Input* pressure_input=inputs->GetInput(PressureEnum); _assert_(pressure_input);
    174176        Input* vx_input=inputs->GetInput(VxEnum);             _assert_(vx_input);
    … …  
    198200        this->inputs->AddInput(new TriaInput(StressTensoryzEnum,&sigma_yz[0],P1Enum));
    199201        this->inputs->AddInput(new TriaInput(StressTensorzzEnum,&sigma_zz[0],P1Enum));
     202
     203        /*Clean up and return*/
     204        delete gauss;
     205}
     206/*}}}*/
     207/*FUNCTION Tria::ComputeDeviatoricStressTensor {{{*/
     208void  Tria::ComputeDeviatoricStressTensor(){
     209
     210        IssmDouble  xyz_list[NUMVERTICES][3];
     211        IssmDouble  viscosity;
     212        IssmDouble  epsilon[3]; /* epsilon=[exx,eyy,exy];*/
     213        IssmDouble  tau_xx[NUMVERTICES];
     214        IssmDouble      tau_yy[NUMVERTICES];
     215        IssmDouble      tau_zz[NUMVERTICES]={0,0,0};
     216        IssmDouble  tau_xy[NUMVERTICES];
     217        IssmDouble      tau_xz[NUMVERTICES]={0,0,0};
     218        IssmDouble      tau_yz[NUMVERTICES]={0,0,0};
     219        GaussTria*  gauss=NULL;
     220        int meshtype,dim=2;
     221
     222        /* Get node coordinates and dof list: */
     223        ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
     224
     225        /*Retrieve all inputs we will be needing: */
     226        this->FindParam(&meshtype,MeshTypeEnum);
     227        if(meshtype!=Mesh2DhorizontalEnum) _error_("deviatoric stress tensor calculation not implemented for mesh of type " <<EnumToStringx(meshtype));
     228        Input* vx_input=inputs->GetInput(VxEnum);             _assert_(vx_input);
     229        Input* vy_input=inputs->GetInput(VyEnum);             _assert_(vy_input);
     230
     231        /* Start looping on the number of vertices: */
     232        gauss=new GaussTria();
     233        for (int iv=0;iv<NUMVERTICES;iv++){
     234                gauss->GaussVertex(iv);
     235
     236                /*Compute strain rate and viscosity: */
     237                this->StrainRateSSA(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input);
     238                this->ViscositySSA(&viscosity,dim,&xyz_list[0][0],gauss,vx_input,vy_input);
     239
     240                /*Compute Stress*/
     241                tau_xx[iv]=2*viscosity*epsilon[0]; // tau = nu eps
     242                tau_yy[iv]=2*viscosity*epsilon[1];
     243                tau_xy[iv]=2*viscosity*epsilon[2];
     244        }
     245
     246        /*Add Stress tensor components into inputs*/
     247        this->inputs->AddInput(new TriaInput(DeviatoricStressxxEnum,&tau_xx[0],P1Enum));
     248        this->inputs->AddInput(new TriaInput(DeviatoricStressxyEnum,&tau_xy[0],P1Enum));
     249        this->inputs->AddInput(new TriaInput(DeviatoricStressxzEnum,&tau_xz[0],P1Enum));
     250        this->inputs->AddInput(new TriaInput(DeviatoricStressyyEnum,&tau_yy[0],P1Enum));
     251        this->inputs->AddInput(new TriaInput(DeviatoricStressyzEnum,&tau_yz[0],P1Enum));
     252        this->inputs->AddInput(new TriaInput(DeviatoricStresszzEnum,&tau_zz[0],P1Enum));
    200253
    201254        /*Clean up and return*/

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

    r17671 r17680  
    5454                void        ComputeStrainRate(Vector<IssmDouble>* eps);
    5555                void        ComputeStressTensor();
     56                void        ComputeDeviatoricStressTensor();
    5657                void        ComputeSurfaceNormalVelocity();
    5758                void        Configure(Elements* elements,Loads* loads,Nodes* nodesin,Vertices* verticesin,Materials* materials,Parameters* parameters);
Note: See TracChangeset for help on using the changeset viewer.