Index: /issm/trunk-jpl/src/c/analyses/DamageEvolutionAnalysis.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/DamageEvolutionAnalysis.cpp (revision 16847) +++ /issm/trunk-jpl/src/c/analyses/DamageEvolutionAnalysis.cpp (revision 16848) @@ -99,5 +99,61 @@ }/*}}}*/ ElementVector* DamageEvolutionAnalysis::CreatePVector(Element* element){/*{{{*/ -_error_("not implemented yet"); + + /*Intermediaries*/ + int meshtype; + Element* basalelement; + + /*Get basal element*/ + element->FindParam(&meshtype,MeshTypeEnum); + switch(meshtype){ + case Mesh2DhorizontalEnum: + basalelement = element; + break; + case Mesh3DEnum: + if(!element->IsOnBed()) return NULL; + basalelement = element->SpawnBasalElement(); + break; + default: _error_("mesh "<GetNumberOfNodes(); + + /*Initialize Element vector and other vectors*/ + ElementVector* pe = element->NewElementVector(); + IssmDouble* basis = xNew(numnodes); + + /*Retrieve all inputs and parameters*/ + element->GetVerticesCoordinates(&xyz_list); + element->FindParam(&dt,TimesteppingTimeStepEnum); + this->CreateDamageFInput(element); + Input* damaged_input = element->GetMaterialInput(DamageDbarEnum); _assert_(damaged_input); + Input* damagef_input = element->GetMaterialInput(DamageFEnum); _assert_(damagef_input); + + /* Start looping on the number of gaussian points: */ + Gauss* gauss=element->NewGauss(2); + for(int ig=gauss->begin();igend();ig++){ + gauss->GaussPoint(ig); + + element->JacobianDeterminant(&Jdet,xyz_list,gauss); + element->NodalFunctions(basis,gauss); + + damaged_input->GetInputValue(&damage,gauss); + damagef_input->GetInputValue(&f,gauss); + + for(int i=0;ivalues[i]+=Jdet*gauss->weight*(damage+dt*f)*basis[i]; + } + + /*Clean up and return*/ + xDelete(xyz_list); + xDelete(basis); + if(meshtype!=Mesh2DhorizontalEnum){basalelement->DeleteMaterials(); delete basalelement;}; + delete gauss; + return pe; }/*}}}*/ void DamageEvolutionAnalysis::GetSolutionFromInputs(Vector* solution,Element* element){/*{{{*/ @@ -136,2 +192,61 @@ xDelete(doflist); }/*}}}*/ + +/*Intermediaries*/ +void DamageEvolutionAnalysis::CreateDamageFInput(Element* element){/*{{{*/ + + /*Intermediaries */ + IssmDouble c1,c2,c3,healing,stress_threshold; + IssmDouble s_xx,s_xy,s_yy; + IssmDouble J2s,Xis,Psi,PosPsi,NegPsi; + IssmDouble damage,sigma_xx,sigma_xy,sigma_yy; + + /*Fetch number of vertices and allocate output*/ + int numvertices = element->GetNumberOfVertices(); + IssmDouble* f = xNew(numvertices); + + /*retrieve parameters:*/ + element->FindParam(&c1,DamageC1Enum); + element->FindParam(&c2,DamageC2Enum); + element->FindParam(&c3,DamageC3Enum); + element->FindParam(&healing,DamageHealingEnum); + element->FindParam(&stress_threshold,DamageStressThresholdEnum); + + /*Compute stress tensor: */ + element->ComputeStressTensor(); + + /*retrieve what we need: */ + Input* sigma_xx_input = element->GetInput(StressTensorxxEnum); _assert_(sigma_xx_input); + Input* sigma_xy_input = element->GetInput(StressTensorxyEnum); _assert_(sigma_xy_input); + Input* sigma_yy_input = element->GetInput(StressTensoryyEnum); _assert_(sigma_yy_input); + Input* damage_input = element->GetMaterialInput(DamageDbarEnum); _assert_(damage_input); + + /*Damage evolution z mapping: */ + Gauss* gauss=element->NewGauss(); + for (int iv=0;ivGaussVertex(iv); + + damage_input->GetInputValue(&damage,gauss); + sigma_xx_input->GetInputValue(&sigma_xx,gauss); + sigma_xy_input->GetInputValue(&sigma_xy,gauss); + sigma_yy_input->GetInputValue(&sigma_yy,gauss); + + s_xx=sigma_xx/(1.-damage); + s_xy=sigma_xy/(1.-damage); + s_yy=sigma_yy/(1.-damage); + + J2s=1./sqrt(2.)*sqrt(s_xx*s_xx + s_yy*s_yy + s_xy*s_xy); + Xis=sqrt(3.0)*J2s; + Psi=Xis-stress_threshold; + PosPsi=max(Psi,0.); + NegPsi=max(-Psi,0.); + + f[iv]= c1*(pow(PosPsi,c2) - healing*pow(NegPsi,c2))*pow((1.-damage),-c3); + } + + /*Add input*/ + element->AddMaterialInput(DamageFEnum,f,P1Enum); + + /*Clean up and return*/ + delete gauss; +}/*}}}*/ Index: /issm/trunk-jpl/src/c/analyses/DamageEvolutionAnalysis.h =================================================================== --- /issm/trunk-jpl/src/c/analyses/DamageEvolutionAnalysis.h (revision 16847) +++ /issm/trunk-jpl/src/c/analyses/DamageEvolutionAnalysis.h (revision 16848) @@ -25,4 +25,7 @@ void GetSolutionFromInputs(Vector* solution,Element* element); void InputUpdateFromSolution(IssmDouble* solution,Element* element); + + /*Intermediaries*/ + void CreateDamageFInput(Element* element); }; #endif Index: /issm/trunk-jpl/src/c/classes/Elements/Element.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Element.h (revision 16847) +++ /issm/trunk-jpl/src/c/classes/Elements/Element.h (revision 16848) @@ -107,4 +107,5 @@ virtual void GetInputListOnVertices(IssmDouble* pvalue,int enumtype,IssmDouble defaultvalue)=0; virtual Input* GetInput(int inputenum)=0; + virtual Input* GetMaterialInput(int inputenum)=0; virtual void GetInputValue(IssmDouble* pvalue,Node* node,int enumtype)=0; virtual void GetInputValue(bool* pvalue,int enum_type)=0; @@ -122,4 +123,5 @@ virtual void ComputeBasalStress(Vector* sigma_b)=0; virtual void ComputeStrainRate(Vector* eps)=0; + virtual void ComputeStressTensor(void)=0; virtual void ResultInterpolation(int* pinterpolation,int output_enum)=0; virtual void ResultToVector(Vector* vector,int output_enum)=0; Index: /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp (revision 16847) +++ /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp (revision 16848) @@ -1312,4 +1312,9 @@ Input* Penta::GetInput(int inputenum){ return inputs->GetInput(inputenum); +} +/*}}}*/ +/*FUNCTION Penta::GetMaterialInput(int inputenum) {{{*/ +Input* Penta::GetMaterialInput(int inputenum){ + return this->material->inputs->GetInput(inputenum); } /*}}}*/ Index: /issm/trunk-jpl/src/c/classes/Elements/Penta.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Penta.h (revision 16847) +++ /issm/trunk-jpl/src/c/classes/Elements/Penta.h (revision 16848) @@ -224,4 +224,5 @@ int GetElementType(void); Input* GetInput(int inputenum); + Input* GetMaterialInput(int inputenum); void GetInputListOnVertices(IssmDouble* pvalue,int enumtype); void GetInputListOnVertices(IssmDouble* pvalue,int enumtype,IssmDouble defaultvalue); Index: /issm/trunk-jpl/src/c/classes/Elements/Seg.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Seg.h (revision 16847) +++ /issm/trunk-jpl/src/c/classes/Elements/Seg.h (revision 16848) @@ -137,4 +137,5 @@ int VelocityInterpolation(void){_error_("not implemented yet");}; Input* GetInput(int inputenum){_error_("not implemented yet");}; + Input* GetMaterialInput(int inputenum){_error_("not implemented yet");}; void GetInputListOnVertices(IssmDouble* pvalue,int enumtype){_error_("not implemented yet");}; void GetInputListOnVertices(IssmDouble* pvalue,int enumtype,IssmDouble defaultvalue){_error_("not implemented yet");}; Index: /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp (revision 16847) +++ /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp (revision 16848) @@ -1342,4 +1342,9 @@ Input* Tria::GetInput(int inputenum){ return inputs->GetInput(inputenum); +} +/*}}}*/ +/*FUNCTION Tria::GetMaterialInput(int inputenum) {{{*/ +Input* Tria::GetMaterialInput(int inputenum){ + return this->material->inputs->GetInput(inputenum); } /*}}}*/ Index: /issm/trunk-jpl/src/c/classes/Elements/Tria.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Tria.h (revision 16847) +++ /issm/trunk-jpl/src/c/classes/Elements/Tria.h (revision 16848) @@ -265,4 +265,5 @@ IssmDouble GetMaterialParameter(int enum_in); Input* GetInput(int inputenum); + Input* GetMaterialInput(int inputenum); void GetInputListOnVertices(IssmDouble* pvalue,int enumtype); void GetInputListOnVertices(IssmDouble* pvalue,int enumtype,IssmDouble defaultvalue); Index: /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h (revision 16847) +++ /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h (revision 16848) @@ -161,4 +161,5 @@ MaterialsRheologyNEnum, DamageDEnum, + DamageFEnum, QmuDamageDEnum, DamageDbarEnum, Index: /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp (revision 16847) +++ /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp (revision 16848) @@ -169,4 +169,5 @@ case MaterialsRheologyNEnum : return "MaterialsRheologyN"; case DamageDEnum : return "DamageD"; + case DamageFEnum : return "DamageF"; case QmuDamageDEnum : return "QmuDamageD"; case DamageDbarEnum : return "DamageDbar"; Index: /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp (revision 16847) +++ /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp (revision 16848) @@ -172,4 +172,5 @@ else if (strcmp(name,"MaterialsRheologyN")==0) return MaterialsRheologyNEnum; else if (strcmp(name,"DamageD")==0) return DamageDEnum; + else if (strcmp(name,"DamageF")==0) return DamageFEnum; else if (strcmp(name,"QmuDamageD")==0) return QmuDamageDEnum; else if (strcmp(name,"DamageDbar")==0) return DamageDbarEnum; @@ -259,9 +260,9 @@ else if (strcmp(name,"SteadystateMaxiter")==0) return SteadystateMaxiterEnum; else if (strcmp(name,"SteadystateNumRequestedOutputs")==0) return SteadystateNumRequestedOutputsEnum; - else if (strcmp(name,"SteadystateReltol")==0) return SteadystateReltolEnum; else stage=3; } if(stage==3){ - if (strcmp(name,"SteadystateRequestedOutputs")==0) return SteadystateRequestedOutputsEnum; + if (strcmp(name,"SteadystateReltol")==0) return SteadystateReltolEnum; + else if (strcmp(name,"SteadystateRequestedOutputs")==0) return SteadystateRequestedOutputsEnum; else if (strcmp(name,"Surface")==0) return SurfaceEnum; else if (strcmp(name,"SurfaceforcingsPrecipitation")==0) return SurfaceforcingsPrecipitationEnum; @@ -382,9 +383,9 @@ else if (strcmp(name,"DoubleMatArrayParam")==0) return DoubleMatArrayParamEnum; else if (strcmp(name,"DoubleMatParam")==0) return DoubleMatParamEnum; - else if (strcmp(name,"DoubleParam")==0) return DoubleParamEnum; else stage=4; } if(stage==4){ - if (strcmp(name,"DoubleVecParam")==0) return DoubleVecParamEnum; + if (strcmp(name,"DoubleParam")==0) return DoubleParamEnum; + else if (strcmp(name,"DoubleVecParam")==0) return DoubleVecParamEnum; else if (strcmp(name,"Element")==0) return ElementEnum; else if (strcmp(name,"ExternalResult")==0) return ExternalResultEnum; @@ -505,9 +506,9 @@ else if (strcmp(name,"Waterfraction")==0) return WaterfractionEnum; else if (strcmp(name,"Watercolumn")==0) return WatercolumnEnum; - else if (strcmp(name,"BasalFriction")==0) return BasalFrictionEnum; else stage=5; } if(stage==5){ - if (strcmp(name,"ViscousHeating")==0) return ViscousHeatingEnum; + if (strcmp(name,"BasalFriction")==0) return BasalFrictionEnum; + else if (strcmp(name,"ViscousHeating")==0) return ViscousHeatingEnum; else if (strcmp(name,"QmuTemperature")==0) return QmuTemperatureEnum; else if (strcmp(name,"HydrologyWaterVx")==0) return HydrologyWaterVxEnum; @@ -628,9 +629,9 @@ else if (strcmp(name,"Gsl")==0) return GslEnum; else if (strcmp(name,"Option")==0) return OptionEnum; - else if (strcmp(name,"GenericOption")==0) return GenericOptionEnum; else stage=6; } if(stage==6){ - if (strcmp(name,"OptionCell")==0) return OptionCellEnum; + if (strcmp(name,"GenericOption")==0) return GenericOptionEnum; + else if (strcmp(name,"OptionCell")==0) return OptionCellEnum; else if (strcmp(name,"OptionStruct")==0) return OptionStructEnum; else if (strcmp(name,"Paterson")==0) return PatersonEnum; Index: /issm/trunk-jpl/src/m/enum/DamageFEnum.m =================================================================== --- /issm/trunk-jpl/src/m/enum/DamageFEnum.m (revision 16848) +++ /issm/trunk-jpl/src/m/enum/DamageFEnum.m (revision 16848) @@ -0,0 +1,11 @@ +function macro=DamageFEnum() +%DAMAGEFENUM - Enum of DamageF +% +% WARNING: DO NOT MODIFY THIS FILE +% this file has been automatically generated by src/c/shared/Enum/Synchronize.sh +% Please read src/c/shared/Enum/README for more information +% +% Usage: +% macro=DamageFEnum() + +macro=StringToEnum('DamageF'); Index: /issm/trunk-jpl/src/m/enum/EnumDefinitions.py =================================================================== --- /issm/trunk-jpl/src/m/enum/EnumDefinitions.py (revision 16847) +++ /issm/trunk-jpl/src/m/enum/EnumDefinitions.py (revision 16848) @@ -161,4 +161,5 @@ def MaterialsRheologyNEnum(): return StringToEnum("MaterialsRheologyN")[0] def DamageDEnum(): return StringToEnum("DamageD")[0] +def DamageFEnum(): return StringToEnum("DamageF")[0] def QmuDamageDEnum(): return StringToEnum("QmuDamageD")[0] def DamageDbarEnum(): return StringToEnum("DamageDbar")[0]