Index: /issm/trunk-jpl/src/c/classes/Elements/Element.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Element.cpp (revision 18193) +++ /issm/trunk-jpl/src/c/classes/Elements/Element.cpp (revision 18194) @@ -41,4 +41,129 @@ this->inputs->AddInput(input_in); }/*}}}*/ +void Element::ComputeNewDamage(){/*{{{*/ + + IssmDouble *xyz_list=NULL; + IssmDouble eps_xx,eps_xy,eps_yy,eps_xz,eps_yz,eps_zz,eps_eff; + IssmDouble epsmin=1.e-27; + int dim; + + /*Retrieve all inputs we will be needing: */ + /* TODO: retrieve parameters eps_0 and eps_f and input DamageD(bar?) */ + this->GetVerticesCoordinates(&xyz_list); + this->ComputeStrainRate(); + parameters->FindParam(&dim,DomainDimensionEnum); + Input* eps_xx_input=this->GetInput(StrainRatexxEnum); _assert_(eps_xx_input); + Input* eps_yy_input=this->GetInput(StrainRateyyEnum); _assert_(eps_yy_input); + Input* eps_xy_input=this->GetInput(StrainRatexyEnum); _assert_(eps_xy_input); + Input* eps_xz_input=NULL; + Input* eps_yz_input=NULL; + Input* eps_zz_input=NULL; + if(dim==3){ + eps_xz_input=this->GetInput(StrainRatexzEnum); _assert_(eps_xz_input); + eps_yz_input=this->GetInput(StrainRateyzEnum); _assert_(eps_yz_input); + eps_zz_input=this->GetInput(StrainRatezzEnum); _assert_(eps_zz_input); + } + + /* Start looping on the number of vertices: */ + Gauss* gauss=this->NewGauss(); + int numvertices = this->GetNumberOfVertices(); + for (int iv=0;ivGaussVertex(iv); + + eps_xx_input->GetInputValue(&eps_xx,gauss); + eps_yy_input->GetInputValue(&eps_yy,gauss); + eps_xy_input->GetInputValue(&eps_xy,gauss); + if(dim==3){ + eps_xz_input->GetInputValue(&eps_xz,gauss); + eps_yz_input->GetInputValue(&eps_yz,gauss); + eps_zz_input->GetInputValue(&eps_zz,gauss); + } + else{eps_xz=0; eps_yz=0; eps_zz=0;} + + /* eps_eff^2 = exx^2 + eyy^2 + exy^2 + exz^2 + eyz^2 + exx*eyy */ + eps_eff=sqrt(eps_xx*eps_xx+eps_yy*eps_yy+eps_xy*eps_xy+eps_xz*eps_xz+eps_yz*eps_yz+eps_xx*eps_xx*epsmin*epsmin); + + /*TODO: compute kappa from initial D, then compute new D */ + + } + + /* TODO: add newdamage input to DamageEnum and NewDamageEnum */ + + /*Clean up and return*/ + xDelete(xyz_list); + delete gauss; + +}/*}}}*/ +void Element::ComputeStrainRate(){/*{{{*/ + + int dim; + IssmDouble *xyz_list = NULL; + IssmDouble epsilon[6]; + + /*Retrieve all inputs we will be needing: */ + this->GetVerticesCoordinates(&xyz_list); + parameters->FindParam(&dim,DomainDimensionEnum); + Input* vx_input=this->GetInput(VxEnum); _assert_(vx_input); + Input* vy_input=this->GetInput(VyEnum); _assert_(vy_input); + Input* vz_input=NULL; + if(dim==3){vz_input=this->GetInput(VzEnum); _assert_(vz_input);} + + /*Allocate arrays*/ + int numvertices = this->GetNumberOfVertices(); + IssmDouble* eps_xx = xNew(numvertices); + IssmDouble* eps_yy = xNew(numvertices); + IssmDouble* eps_zz = xNew(numvertices); + IssmDouble* eps_xy = xNew(numvertices); + IssmDouble* eps_xz = xNew(numvertices); + IssmDouble* eps_yz = xNew(numvertices); + + /* Start looping on the number of vertices: */ + Gauss* gauss=this->NewGauss(); + for (int iv=0;ivGaussVertex(iv); + + /*Compute strain rate viscosity and pressure: */ + if(dim==2) + this->StrainRateSSA(&epsilon[0],xyz_list,gauss,vx_input,vy_input); + else + this->StrainRateFS(&epsilon[0],xyz_list,gauss,vx_input,vy_input,vz_input); + + if(dim==2){ + /* epsilon=[exx,eyy,exy];*/ + eps_xx[iv]=epsilon[0]; + eps_yy[iv]=epsilon[1]; + eps_xy[iv]=epsilon[2]; + } + else{ + /*epsilon=[exx eyy ezz exy exz eyz]*/ + eps_xx[iv]=epsilon[0]; + eps_yy[iv]=epsilon[1]; + eps_zz[iv]=epsilon[2]; + eps_xy[iv]=epsilon[3]; + eps_xz[iv]=epsilon[4]; + eps_yz[iv]=epsilon[5]; + } + } + + /*Add Stress tensor components into inputs*/ + this->AddInput(StrainRatexxEnum,eps_xx,P1Enum); + this->AddInput(StrainRatexyEnum,eps_xy,P1Enum); + this->AddInput(StrainRatexzEnum,eps_xz,P1Enum); + this->AddInput(StrainRateyyEnum,eps_yy,P1Enum); + this->AddInput(StrainRateyzEnum,eps_yz,P1Enum); + this->AddInput(StrainRatezzEnum,eps_zz,P1Enum); + + /*Clean up and return*/ + delete gauss; + xDelete(xyz_list); + xDelete(eps_xx); + xDelete(eps_yy); + xDelete(eps_zz); + xDelete(eps_xy); + xDelete(eps_xz); + xDelete(eps_yz); + +} +/*}}}*/ void Element::CoordinateSystemTransform(IssmDouble** ptransform,Node** nodes_list,int numnodes,int* cs_array){/*{{{*/ @@ -113,75 +238,4 @@ /*Assign output pointer*/ *ptransform=transform; -} -/*}}}*/ -void Element::ComputeStrainRate(){/*{{{*/ - - int dim; - IssmDouble *xyz_list = NULL; - IssmDouble epsilon[6]; - - /*Retrieve all inputs we will be needing: */ - this->GetVerticesCoordinates(&xyz_list); - parameters->FindParam(&dim,DomainDimensionEnum); - Input* vx_input=this->GetInput(VxEnum); _assert_(vx_input); - Input* vy_input=this->GetInput(VyEnum); _assert_(vy_input); - Input* vz_input=NULL; - if(dim==3){vz_input=this->GetInput(VzEnum); _assert_(vz_input);} - - /*Allocate arrays*/ - int numvertices = this->GetNumberOfVertices(); - IssmDouble* eps_xx = xNew(numvertices); - IssmDouble* eps_yy = xNew(numvertices); - IssmDouble* eps_zz = xNew(numvertices); - IssmDouble* eps_xy = xNew(numvertices); - IssmDouble* eps_xz = xNew(numvertices); - IssmDouble* eps_yz = xNew(numvertices); - - /* Start looping on the number of vertices: */ - Gauss* gauss=this->NewGauss(); - for (int iv=0;ivGaussVertex(iv); - - /*Compute strain rate viscosity and pressure: */ - if(dim==2) - this->StrainRateSSA(&epsilon[0],xyz_list,gauss,vx_input,vy_input); - else - this->StrainRateFS(&epsilon[0],xyz_list,gauss,vx_input,vy_input,vz_input); - - if(dim==2){ - /* epsilon=[exx,eyy,exy];*/ - eps_xx[iv]=epsilon[0]; - eps_yy[iv]=epsilon[1]; - eps_xy[iv]=epsilon[2]; - } - else{ - /*epsilon=[exx eyy ezz exy exz eyz]*/ - eps_xx[iv]=epsilon[0]; - eps_yy[iv]=epsilon[1]; - eps_zz[iv]=epsilon[2]; - eps_xy[iv]=epsilon[3]; - eps_xz[iv]=epsilon[4]; - eps_yz[iv]=epsilon[5]; - } - } - - /*Add Stress tensor components into inputs*/ - this->AddInput(StrainRatexxEnum,eps_xx,P1Enum); - this->AddInput(StrainRatexyEnum,eps_xy,P1Enum); - this->AddInput(StrainRatexzEnum,eps_xz,P1Enum); - this->AddInput(StrainRateyyEnum,eps_yy,P1Enum); - this->AddInput(StrainRateyzEnum,eps_yz,P1Enum); - this->AddInput(StrainRatezzEnum,eps_zz,P1Enum); - - /*Clean up and return*/ - delete gauss; - xDelete(xyz_list); - xDelete(eps_xx); - xDelete(eps_yy); - xDelete(eps_zz); - xDelete(eps_xy); - xDelete(eps_xz); - xDelete(eps_yz); - } /*}}}*/ @@ -1053,4 +1107,8 @@ input=this->inputs->GetInput(output_enum); break; + case NewDamageEnum: + this->ComputeNewDamage(); + input=this->inputs->GetInput(output_enum); + break; default: _error_("input "<