Changeset 16848


Ignore:
Timestamp:
11/21/13 09:21:32 (11 years ago)
Author:
Mathieu Morlighem
Message:

CHG: done with Damage evolution

Location:
issm/trunk-jpl/src
Files:
1 added
12 edited

Legend:

Unmodified
Added
Removed

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

    r16782 r16848  
    9999}/*}}}*/
    100100ElementVector* DamageEvolutionAnalysis::CreatePVector(Element* element){/*{{{*/
    101 _error_("not implemented yet");
     101
     102        /*Intermediaries*/
     103        int      meshtype;
     104        Element* basalelement;
     105
     106        /*Get basal element*/
     107        element->FindParam(&meshtype,MeshTypeEnum);
     108        switch(meshtype){
     109                case Mesh2DhorizontalEnum:
     110                        basalelement = element;
     111                        break;
     112                case Mesh3DEnum:
     113                        if(!element->IsOnBed()) return NULL;
     114                        basalelement = element->SpawnBasalElement();
     115                        break;
     116                default: _error_("mesh "<<EnumToStringx(meshtype)<<" not supported yet");
     117        }
     118
     119        /*Intermediaries */
     120        IssmDouble  Jdet,dt;
     121        IssmDouble  f,damage;
     122        IssmDouble* xyz_list = NULL;
     123
     124        /*Fetch number of nodes and dof for this finite element*/
     125        int numnodes = element->GetNumberOfNodes();
     126
     127        /*Initialize Element vector and other vectors*/
     128        ElementVector* pe    = element->NewElementVector();
     129        IssmDouble*    basis = xNew<IssmDouble>(numnodes);
     130
     131        /*Retrieve all inputs and parameters*/
     132        element->GetVerticesCoordinates(&xyz_list);
     133        element->FindParam(&dt,TimesteppingTimeStepEnum);
     134        this->CreateDamageFInput(element);
     135        Input* damaged_input = element->GetMaterialInput(DamageDbarEnum); _assert_(damaged_input);
     136        Input* damagef_input = element->GetMaterialInput(DamageFEnum);    _assert_(damagef_input);
     137
     138        /* Start  looping on the number of gaussian points: */
     139        Gauss* gauss=element->NewGauss(2);
     140        for(int ig=gauss->begin();ig<gauss->end();ig++){
     141                gauss->GaussPoint(ig);
     142
     143                element->JacobianDeterminant(&Jdet,xyz_list,gauss);
     144                element->NodalFunctions(basis,gauss);
     145
     146                damaged_input->GetInputValue(&damage,gauss);
     147                damagef_input->GetInputValue(&f,gauss);
     148
     149                for(int i=0;i<numnodes;i++) pe->values[i]+=Jdet*gauss->weight*(damage+dt*f)*basis[i];
     150        }
     151
     152        /*Clean up and return*/
     153        xDelete<IssmDouble>(xyz_list);
     154        xDelete<IssmDouble>(basis);
     155        if(meshtype!=Mesh2DhorizontalEnum){basalelement->DeleteMaterials(); delete basalelement;};
     156        delete gauss;
     157        return pe;
    102158}/*}}}*/
    103159void DamageEvolutionAnalysis::GetSolutionFromInputs(Vector<IssmDouble>* solution,Element* element){/*{{{*/
    … …  
    136192        xDelete<int>(doflist);
    137193}/*}}}*/
     194
     195/*Intermediaries*/
     196void DamageEvolutionAnalysis::CreateDamageFInput(Element* element){/*{{{*/
     197
     198        /*Intermediaries */
     199        IssmDouble c1,c2,c3,healing,stress_threshold;
     200        IssmDouble s_xx,s_xy,s_yy;
     201        IssmDouble J2s,Xis,Psi,PosPsi,NegPsi;
     202        IssmDouble damage,sigma_xx,sigma_xy,sigma_yy;
     203
     204        /*Fetch number of vertices and allocate output*/
     205        int numvertices = element->GetNumberOfVertices();
     206        IssmDouble* f   = xNew<IssmDouble>(numvertices);
     207
     208        /*retrieve parameters:*/
     209        element->FindParam(&c1,DamageC1Enum);
     210        element->FindParam(&c2,DamageC2Enum);
     211        element->FindParam(&c3,DamageC3Enum);
     212        element->FindParam(&healing,DamageHealingEnum);
     213        element->FindParam(&stress_threshold,DamageStressThresholdEnum);
     214
     215        /*Compute stress tensor: */
     216        element->ComputeStressTensor();
     217
     218        /*retrieve what we need: */
     219        Input* sigma_xx_input  = element->GetInput(StressTensorxxEnum);     _assert_(sigma_xx_input);
     220        Input* sigma_xy_input  = element->GetInput(StressTensorxyEnum);     _assert_(sigma_xy_input);
     221        Input* sigma_yy_input  = element->GetInput(StressTensoryyEnum);     _assert_(sigma_yy_input);
     222        Input* damage_input    = element->GetMaterialInput(DamageDbarEnum); _assert_(damage_input);
     223
     224        /*Damage evolution z mapping: */
     225        Gauss* gauss=element->NewGauss();
     226        for (int iv=0;iv<numvertices;iv++){
     227                gauss->GaussVertex(iv);
     228               
     229                damage_input->GetInputValue(&damage,gauss);
     230                sigma_xx_input->GetInputValue(&sigma_xx,gauss);
     231                sigma_xy_input->GetInputValue(&sigma_xy,gauss);
     232                sigma_yy_input->GetInputValue(&sigma_yy,gauss);
     233
     234                s_xx=sigma_xx/(1.-damage);
     235                s_xy=sigma_xy/(1.-damage);
     236                s_yy=sigma_yy/(1.-damage);
     237
     238                J2s=1./sqrt(2.)*sqrt(s_xx*s_xx + s_yy*s_yy + s_xy*s_xy);
     239                Xis=sqrt(3.0)*J2s;
     240                Psi=Xis-stress_threshold;
     241                PosPsi=max(Psi,0.);
     242                NegPsi=max(-Psi,0.);
     243
     244                f[iv]= c1*(pow(PosPsi,c2) - healing*pow(NegPsi,c2))*pow((1.-damage),-c3);
     245        }
     246
     247        /*Add input*/
     248        element->AddMaterialInput(DamageFEnum,f,P1Enum);
     249       
     250        /*Clean up and return*/
     251        delete gauss;
     252}/*}}}*/

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

    r16782 r16848  
    2525                void GetSolutionFromInputs(Vector<IssmDouble>* solution,Element* element);
    2626                void InputUpdateFromSolution(IssmDouble* solution,Element* element);
     27
     28                /*Intermediaries*/
     29                void CreateDamageFInput(Element* element);
    2730};
    2831#endif

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

    r16847 r16848  
    107107                virtual void   GetInputListOnVertices(IssmDouble* pvalue,int enumtype,IssmDouble defaultvalue)=0;
    108108                virtual Input* GetInput(int inputenum)=0;
     109                virtual Input* GetMaterialInput(int inputenum)=0;
    109110                virtual void   GetInputValue(IssmDouble* pvalue,Node* node,int enumtype)=0;
    110111                virtual void   GetInputValue(bool* pvalue,int enum_type)=0;
    … …  
    122123                virtual void   ComputeBasalStress(Vector<IssmDouble>* sigma_b)=0;
    123124                virtual void   ComputeStrainRate(Vector<IssmDouble>* eps)=0;
     125                virtual void   ComputeStressTensor(void)=0;
    124126                virtual void   ResultInterpolation(int* pinterpolation,int output_enum)=0;
    125127                virtual void   ResultToVector(Vector<IssmDouble>* vector,int output_enum)=0;

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

    r16847 r16848  
    13121312Input* Penta::GetInput(int inputenum){
    13131313        return inputs->GetInput(inputenum);
     1314}
     1315/*}}}*/
     1316/*FUNCTION Penta::GetMaterialInput(int inputenum) {{{*/
     1317Input* Penta::GetMaterialInput(int inputenum){
     1318        return this->material->inputs->GetInput(inputenum);
    13141319}
    13151320/*}}}*/

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

    r16847 r16848  
    224224                int            GetElementType(void);
    225225                Input*         GetInput(int inputenum);
     226                Input*         GetMaterialInput(int inputenum);
    226227                void           GetInputListOnVertices(IssmDouble* pvalue,int enumtype);
    227228                void           GetInputListOnVertices(IssmDouble* pvalue,int enumtype,IssmDouble defaultvalue);

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

    r16847 r16848  
    137137                int         VelocityInterpolation(void){_error_("not implemented yet");};
    138138                Input*      GetInput(int inputenum){_error_("not implemented yet");};
     139                Input*      GetMaterialInput(int inputenum){_error_("not implemented yet");};
    139140                void        GetInputListOnVertices(IssmDouble* pvalue,int enumtype){_error_("not implemented yet");};
    140141                void        GetInputListOnVertices(IssmDouble* pvalue,int enumtype,IssmDouble defaultvalue){_error_("not implemented yet");};

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

    r16847 r16848  
    13421342Input* Tria::GetInput(int inputenum){
    13431343        return inputs->GetInput(inputenum);
     1344}
     1345/*}}}*/
     1346/*FUNCTION Tria::GetMaterialInput(int inputenum) {{{*/
     1347Input* Tria::GetMaterialInput(int inputenum){
     1348        return this->material->inputs->GetInput(inputenum);
    13441349}
    13451350/*}}}*/

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

    r16847 r16848  
    265265                IssmDouble     GetMaterialParameter(int enum_in);
    266266                Input*         GetInput(int inputenum);
     267                Input*         GetMaterialInput(int inputenum);
    267268                void           GetInputListOnVertices(IssmDouble* pvalue,int enumtype);
    268269                void           GetInputListOnVertices(IssmDouble* pvalue,int enumtype,IssmDouble defaultvalue);

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

    r16787 r16848  
    161161        MaterialsRheologyNEnum,
    162162        DamageDEnum,
     163        DamageFEnum,
    163164        QmuDamageDEnum,
    164165        DamageDbarEnum,

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

    r16787 r16848  
    169169                case MaterialsRheologyNEnum : return "MaterialsRheologyN";
    170170                case DamageDEnum : return "DamageD";
     171                case DamageFEnum : return "DamageF";
    171172                case QmuDamageDEnum : return "QmuDamageD";
    172173                case DamageDbarEnum : return "DamageDbar";

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

    r16787 r16848  
    172172              else if (strcmp(name,"MaterialsRheologyN")==0) return MaterialsRheologyNEnum;
    173173              else if (strcmp(name,"DamageD")==0) return DamageDEnum;
     174              else if (strcmp(name,"DamageF")==0) return DamageFEnum;
    174175              else if (strcmp(name,"QmuDamageD")==0) return QmuDamageDEnum;
    175176              else if (strcmp(name,"DamageDbar")==0) return DamageDbarEnum;
    … …  
    259260              else if (strcmp(name,"SteadystateMaxiter")==0) return SteadystateMaxiterEnum;
    260261              else if (strcmp(name,"SteadystateNumRequestedOutputs")==0) return SteadystateNumRequestedOutputsEnum;
    261               else if (strcmp(name,"SteadystateReltol")==0) return SteadystateReltolEnum;
    262262         else stage=3;
    263263   }
    264264   if(stage==3){
    265               if (strcmp(name,"SteadystateRequestedOutputs")==0) return SteadystateRequestedOutputsEnum;
     265              if (strcmp(name,"SteadystateReltol")==0) return SteadystateReltolEnum;
     266              else if (strcmp(name,"SteadystateRequestedOutputs")==0) return SteadystateRequestedOutputsEnum;
    266267              else if (strcmp(name,"Surface")==0) return SurfaceEnum;
    267268              else if (strcmp(name,"SurfaceforcingsPrecipitation")==0) return SurfaceforcingsPrecipitationEnum;
    … …  
    382383              else if (strcmp(name,"DoubleMatArrayParam")==0) return DoubleMatArrayParamEnum;
    383384              else if (strcmp(name,"DoubleMatParam")==0) return DoubleMatParamEnum;
    384               else if (strcmp(name,"DoubleParam")==0) return DoubleParamEnum;
    385385         else stage=4;
    386386   }
    387387   if(stage==4){
    388               if (strcmp(name,"DoubleVecParam")==0) return DoubleVecParamEnum;
     388              if (strcmp(name,"DoubleParam")==0) return DoubleParamEnum;
     389              else if (strcmp(name,"DoubleVecParam")==0) return DoubleVecParamEnum;
    389390              else if (strcmp(name,"Element")==0) return ElementEnum;
    390391              else if (strcmp(name,"ExternalResult")==0) return ExternalResultEnum;
    … …  
    505506              else if (strcmp(name,"Waterfraction")==0) return WaterfractionEnum;
    506507              else if (strcmp(name,"Watercolumn")==0) return WatercolumnEnum;
    507               else if (strcmp(name,"BasalFriction")==0) return BasalFrictionEnum;
    508508         else stage=5;
    509509   }
    510510   if(stage==5){
    511               if (strcmp(name,"ViscousHeating")==0) return ViscousHeatingEnum;
     511              if (strcmp(name,"BasalFriction")==0) return BasalFrictionEnum;
     512              else if (strcmp(name,"ViscousHeating")==0) return ViscousHeatingEnum;
    512513              else if (strcmp(name,"QmuTemperature")==0) return QmuTemperatureEnum;
    513514              else if (strcmp(name,"HydrologyWaterVx")==0) return HydrologyWaterVxEnum;
    … …  
    628629              else if (strcmp(name,"Gsl")==0) return GslEnum;
    629630              else if (strcmp(name,"Option")==0) return OptionEnum;
    630               else if (strcmp(name,"GenericOption")==0) return GenericOptionEnum;
    631631         else stage=6;
    632632   }
    633633   if(stage==6){
    634               if (strcmp(name,"OptionCell")==0) return OptionCellEnum;
     634              if (strcmp(name,"GenericOption")==0) return GenericOptionEnum;
     635              else if (strcmp(name,"OptionCell")==0) return OptionCellEnum;
    635636              else if (strcmp(name,"OptionStruct")==0) return OptionStructEnum;
    636637              else if (strcmp(name,"Paterson")==0) return PatersonEnum;

  • issm/trunk-jpl/src/m/enum/EnumDefinitions.py

    r16787 r16848  
    161161def MaterialsRheologyNEnum(): return StringToEnum("MaterialsRheologyN")[0]
    162162def DamageDEnum(): return StringToEnum("DamageD")[0]
     163def DamageFEnum(): return StringToEnum("DamageF")[0]
    163164def QmuDamageDEnum(): return StringToEnum("QmuDamageD")[0]
    164165def DamageDbarEnum(): return StringToEnum("DamageDbar")[0]
Note: See TracChangeset for help on using the changeset viewer.