Changeset 18470


Ignore:
Timestamp:
08/26/14 16:45:22 (11 years ago)
Author:
srebuffi
Message:

NEW: Added damage evolution in 3d

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

Legend:

Unmodified
Added
Removed

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

    r18277 r18470  
    7575
    7676        iomodel->Constant(&finiteelement,DamageElementinterpEnum);
    77 
    78         if(iomodel->domaintype!=Domain2DhorizontalEnum) iomodel->FetchData(1,MeshVertexonbaseEnum);
    7977        ::CreateNodes(nodes,iomodel,DamageEvolutionAnalysisEnum,finiteelement);
    80         iomodel->DeleteData(1,MeshVertexonbaseEnum);
    8178}/*}}}*/
    8279void DamageEvolutionAnalysis::CreateConstraints(Constraints* constraints,IoModel* iomodel){/*{{{*/
    … …  
    109106        /* Check if ice in element */
    110107        if(!element->IsIceInElement()) return NULL;
    111 
    112108        /*Intermediaries*/
    113         Element*    basalelement;
    114109        int         domaintype,dim;
    115110        int         stabilization;
    116         IssmDouble  Jdet,dt,D_scalar,h;
    117         IssmDouble  vel,vx,vy,dvxdx,dvydy,dvx[2],dvy[2];
     111        IssmDouble  Jdet,dt,D_scalar,h,hx,hy,hz;
     112        IssmDouble  vel,vx,vy,vz,dvxdx,dvydy,dvzdz,dvx[3],dvy[3],dvz[3];
    118113        IssmDouble *xyz_list  = NULL;
    119114
    120         /*Get problem dimension and basal element*/
     115        /*Get problem dimension*/
    121116        element->FindParam(&domaintype,DomainTypeEnum);
    122117        switch(domaintype){
    123                 case Domain2DhorizontalEnum:
    124                         basalelement = element;
    125                         dim = 2;
    126                         break;
    127                 case Domain3DEnum:
    128                         if(!element->IsOnBase()) return NULL;
    129                         basalelement = element->SpawnBasalElement();
    130                         dim = 2;
    131                         break;
    132                 default: _error_("mesh "<<EnumToStringx(domaintype)<<" not supported yet");
     118                case Domain2DhorizontalEnum: dim = 2; break;
     119                case Domain3DEnum:           dim = 3; break;
     120                default: _error_("Not implemented yet");
    133121        }
    134122
    135123        /*Fetch number of nodes and dof for this finite element*/
    136         int numnodes = basalelement->GetNumberOfNodes();
     124        int numnodes = element->GetNumberOfNodes();
    137125
    138126        /*Initialize Element vector*/
    139         ElementMatrix* Ke     = basalelement->NewElementMatrix();
     127        ElementMatrix* Ke     = element->NewElementMatrix();
    140128        IssmDouble*    basis  = xNew<IssmDouble>(numnodes);
    141129        IssmDouble*    B      = xNew<IssmDouble>(dim*numnodes);
    … …  
    144132
    145133        /*Retrieve all inputs and parameters*/
    146         basalelement->GetVerticesCoordinates(&xyz_list);
    147         basalelement->FindParam(&dt,TimesteppingTimeStepEnum);
    148         basalelement->FindParam(&stabilization,DamageStabilizationEnum);
    149         Input* vxaverage_input=NULL;
    150         Input* vyaverage_input=NULL;
    151         if(domaintype==Domain2DhorizontalEnum){
    152                 vxaverage_input=basalelement->GetInput(VxEnum); _assert_(vxaverage_input);
    153                 vyaverage_input=basalelement->GetInput(VyEnum); _assert_(vyaverage_input);
    154         }
    155         else{
    156                 if(dim==1){
    157                         vxaverage_input=element->GetInput(VxEnum); _assert_(vxaverage_input);
    158                 }
    159                 if(dim==2){
    160                         vxaverage_input=element->GetInput(VxAverageEnum); _assert_(vxaverage_input);
    161                         vyaverage_input=element->GetInput(VyAverageEnum); _assert_(vyaverage_input);
    162                 }
    163         }
    164         h=basalelement->CharacteristicLength();
     134        element->GetVerticesCoordinates(&xyz_list);
     135        element->FindParam(&dt,TimesteppingTimeStepEnum);
     136        element->FindParam(&stabilization,DamageStabilizationEnum);
     137        Input* vx_input = element->GetInput(VxEnum); _assert_(vx_input);
     138        Input* vy_input = element->GetInput(VyEnum); _assert_(vy_input);
     139        Input* vz_input = NULL;
     140        if(dim==3){
     141                vz_input=element->GetInput(VzEnum); _assert_(vz_input);
     142        }
     143
     144        if(dim==2) h=element->CharacteristicLength();
    165145
    166146        /* Start  looping on the number of gaussian points: */
    167         Gauss* gauss=basalelement->NewGauss(2);
     147        Gauss* gauss=element->NewGauss(2);
    168148        for(int ig=gauss->begin();ig<gauss->end();ig++){
    169149                gauss->GaussPoint(ig);
    170150
    171                 basalelement->JacobianDeterminant(&Jdet,xyz_list,gauss);
    172                 basalelement->NodalFunctions(basis,gauss);
     151                element->JacobianDeterminant(&Jdet,xyz_list,gauss);
     152                element->NodalFunctions(basis,gauss);
    173153               
    174                 vxaverage_input->GetInputValue(&vx,gauss);
    175                 vxaverage_input->GetInputDerivativeValue(&dvx[0],xyz_list,gauss);
    176                 if(dim==2){
    177                         vyaverage_input->GetInputValue(&vy,gauss);
    178                         vyaverage_input->GetInputDerivativeValue(&dvy[0],xyz_list,gauss);
     154                vx_input->GetInputValue(&vx,gauss);
     155                vx_input->GetInputDerivativeValue(&dvx[0],xyz_list,gauss);
     156                vy_input->GetInputValue(&vy,gauss);
     157                vy_input->GetInputDerivativeValue(&dvy[0],xyz_list,gauss);
     158
     159                if(dim==3){
     160                        vz_input->GetInputValue(&vz,gauss);
     161                        vz_input->GetInputDerivativeValue(&dvz[0],xyz_list,gauss);
    179162                }
    180163
    … …  
    185168                                        &Ke->values[0],1);
    186169
    187                 GetB(B,basalelement,dim,xyz_list,gauss);
    188                 GetBprime(Bprime,basalelement,dim,xyz_list,gauss);
     170                GetB(B,element,dim,xyz_list,gauss);
     171                GetBprime(Bprime,element,dim,xyz_list,gauss);
    189172
    190173                dvxdx=dvx[0];
    191                 if(dim==2) dvydy=dvy[1];
     174                dvydy=dvy[1];
     175                if(dim==3) dvzdz=dvz[2];
    192176                D_scalar=dt*gauss->weight*Jdet;
    193177
    194178                D[0*dim+0]=D_scalar*dvxdx;
    195                 if(dim==2) D[1*dim+1]=D_scalar*dvydy;
     179                D[1*dim+1]=D_scalar*dvydy;
     180                if(dim==3) D[2*dim+2]=D_scalar*dvzdz;
    196181
    197182                TripleMultiply(B,dim,numnodes,1,
    … …  
    201186
    202187                D[0*dim+0]=D_scalar*vx;
    203                 if(dim==2) D[1*dim+1]=D_scalar*vy;
     188                D[1*dim+1]=D_scalar*vy;
     189                if(dim==3) D[2*dim+2]=D_scalar*vz;
    204190
    205191                TripleMultiply(B,dim,numnodes,1,
    … …  
    209195
    210196                if(stabilization==2){
    211                         if(dim==1){
    212                                 vel=fabs(vx)+1.e-8;
    213                                 D[0]=h/(2.0*vel)*vx*vx;
     197                        if(dim==3){
     198                                vel=sqrt(vx*vx+vy*vy+vz*vz)+1.e-8;
     199                                D[0*dim+0]=h/(2.0*vel)*vx*vx;
     200                                D[1*dim+0]=h/(2.0*vel)*vy*vx;
     201                                D[2*dim+0]=h/(2.0*vel)*vz*vx;
     202                                D[0*dim+1]=h/(2.0*vel)*vx*vy;
     203                                D[1*dim+1]=h/(2.0*vel)*vy*vy;
     204                                D[2*dim+1]=h/(2.0*vel)*vy*vz;
     205                                D[0*dim+2]=h/(2.0*vel)*vx*vz;
     206                                D[1*dim+2]=h/(2.0*vel)*vy*vz;
     207                                D[2*dim+2]=h/(2.0*vel)*vz*vz;
    214208                        }
    215209                        else{
    … …  
    223217                }
    224218                else if(stabilization==1){
    225                         vxaverage_input->GetInputAverage(&vx);
    226                         if(dim==2) vyaverage_input->GetInputAverage(&vy);
    227                         D[0*dim+0]=h/2.0*fabs(vx);
    228                         if(dim==2) D[1*dim+1]=h/2.0*fabs(vy);
     219                        if(dim==2){
     220                                vx_input->GetInputAverage(&vx);
     221                                vy_input->GetInputAverage(&vy);
     222                                D[0*dim+0]=h/2.0*fabs(vx);
     223                                D[1*dim+1]=h/2.0*fabs(vy);
     224                        }
     225                        else if(dim==3){
     226                                element->ElementSizes(&hx,&hy,&hz);
     227                                vel=sqrt(vx*vx + vy*vy + vz*vz)+1.e-14;
     228                                h=sqrt( pow(hx*vx/vel,2) + pow(hy*vy/vel,2) + pow(hz*vz/vel,2));
     229                                D[0*dim+0]=h/(2.*vel)*fabs(vx*vx);  D[0*dim+1]=h/(2.*vel)*fabs(vx*vy); D[0*dim+2]=h/(2.*vel)*fabs(vx*vz);
     230                                D[1*dim+0]=h/(2.*vel)*fabs(vy*vx);  D[1*dim+1]=h/(2.*vel)*fabs(vy*vy); D[1*dim+2]=h/(2.*vel)*fabs(vy*vz);
     231                                D[2*dim+0]=h/(2.*vel)*fabs(vz*vx);  D[2*dim+1]=h/(2.*vel)*fabs(vz*vy); D[2*dim+2]=h/(2.*vel)*fabs(vz*vz);
     232                        }
    229233                }
    230234                if(stabilization==1 || stabilization==2){
    231                         if(dim==1) D[0]=D_scalar*D[0];
    232                         else{
     235                        if(dim==2){
    233236                                D[0*dim+0]=D_scalar*D[0*dim+0];
    234237                                D[1*dim+0]=D_scalar*D[1*dim+0];
    … …  
    236239                                D[1*dim+1]=D_scalar*D[1*dim+1];
    237240                        }
    238 
     241                        else if(dim==3){
     242                                D[0*dim+0]=D_scalar*D[0*dim+0];
     243                                D[1*dim+0]=D_scalar*D[1*dim+0];
     244                                D[2*dim+0]=D_scalar*D[2*dim+0];
     245                                D[0*dim+1]=D_scalar*D[0*dim+1];
     246                                D[1*dim+1]=D_scalar*D[1*dim+1];
     247                                D[2*dim+1]=D_scalar*D[2*dim+1];
     248                                D[0*dim+2]=D_scalar*D[0*dim+2];
     249                                D[1*dim+2]=D_scalar*D[1*dim+2];
     250                                D[2*dim+2]=D_scalar*D[2*dim+2];
     251                        }
    239252                        TripleMultiply(Bprime,dim,numnodes,1,
    240253                                                D,dim,dim,0,
    … …  
    252265        xDelete<IssmDouble>(D);
    253266        delete gauss;
    254         if(domaintype!=Domain2DhorizontalEnum){basalelement->DeleteMaterials(); delete basalelement;};
    255267        return Ke;
    256268}/*}}}*/
    … …  
    262274        /*Intermediaries*/
    263275        int      domaintype,damagelaw;
    264         Element* basalelement;
    265276        IssmDouble  Jdet,dt;
    266277        IssmDouble  f,damage;
    267278        IssmDouble* xyz_list = NULL;
    268 
    269         /*Get basal element*/
     279        /*Get element*/
    270280        element->FindParam(&domaintype,DomainTypeEnum);
    271         switch(domaintype){
    272                 case Domain2DhorizontalEnum:
    273                         basalelement = element;
     281
     282        /*Fetch number of nodes and dof for this finite element*/
     283        int numnodes = element->GetNumberOfNodes();
     284
     285        /*Initialize Element vector and other vectors*/
     286        ElementVector* pe    = element->NewElementVector();
     287        IssmDouble*    basis = xNew<IssmDouble>(numnodes);
     288
     289        /*Retrieve all inputs and parameters*/
     290        element->GetVerticesCoordinates(&xyz_list);
     291        element->FindParam(&dt,TimesteppingTimeStepEnum);
     292        element->FindParam(&damagelaw,DamageLawEnum);
     293        switch(damagelaw){
     294                case 1:
     295                        this->CreateDamageFInputPralong(element);
    274296                        break;
    275                 case Domain3DEnum:
    276                         if(!element->IsOnBase()) return NULL;
    277                         basalelement = element->SpawnBasalElement();
     297                case 2:
     298                        this->CreateDamageFInputPralong(element);
    278299                        break;
    279                 default: _error_("mesh "<<EnumToStringx(domaintype)<<" not supported yet");
    280         }
    281 
    282 
    283         /*Fetch number of nodes and dof for this finite element*/
    284         int numnodes = basalelement->GetNumberOfNodes();
    285 
    286         /*Initialize Element vector and other vectors*/
    287         ElementVector* pe    = basalelement->NewElementVector();
    288         IssmDouble*    basis = xNew<IssmDouble>(numnodes);
    289 
    290         /*Retrieve all inputs and parameters*/
    291         basalelement->GetVerticesCoordinates(&xyz_list);
    292         basalelement->FindParam(&dt,TimesteppingTimeStepEnum);
    293         basalelement->FindParam(&damagelaw,DamageLawEnum);
    294         if(damagelaw==1 | damagelaw==2){
    295                 this->CreateDamageFInputPralong(basalelement);
    296         }
    297         else if(damagelaw==3){
    298                 this->CreateDamageFInputExp(basalelement);
    299         }
    300         Input* damaged_input = basalelement->GetInput(DamageDEnum); _assert_(damaged_input);
    301         Input* damagef_input = basalelement->GetInput(DamageFEnum); _assert_(damagef_input);
     300                case 3:
     301                        this->CreateDamageFInputExp(element);
     302                        break;
     303                default:
     304                        _error_("not implemented yet");
     305        }
     306
     307        Input* damaged_input = NULL;
     308        Input* damagef_input = element->GetInput(DamageFEnum); _assert_(damagef_input);
     309        if(domaintype==Domain2DhorizontalEnum){
     310                damaged_input = element->GetInput(DamageDbarEnum); _assert_(damaged_input);
     311        }
     312        else{
     313                damaged_input = element->GetInput(DamageDEnum); _assert_(damaged_input);
     314        }
    302315
    303316
    304317        /* Start  looping on the number of gaussian points: */
    305         Gauss* gauss=basalelement->NewGauss(2);
     318        Gauss* gauss=element->NewGauss(2);
    306319        for(int ig=gauss->begin();ig<gauss->end();ig++){
    307320                gauss->GaussPoint(ig);
    308321
    309                 basalelement->JacobianDeterminant(&Jdet,xyz_list,gauss);
    310                 basalelement->NodalFunctions(basis,gauss);
     322                element->JacobianDeterminant(&Jdet,xyz_list,gauss);
     323                element->NodalFunctions(basis,gauss);
    311324
    312325                damaged_input->GetInputValue(&damage,gauss);
    … …  
    317330                }
    318331        }
    319 
    320332        /*Clean up and return*/
    321333        xDelete<IssmDouble>(xyz_list);
    322334        xDelete<IssmDouble>(basis);
    323         if(domaintype!=Domain2DhorizontalEnum){basalelement->DeleteMaterials(); delete basalelement;};
    324335        delete gauss;
    325336        return pe;
    … …  
    393404        IssmDouble  max_damage;
    394405        int                     *doflist = NULL;
    395         Element*   basalelement=NULL;
    396406
    397407        element->FindParam(&domaintype,DomainTypeEnum);
    398         if(domaintype!=Domain2DhorizontalEnum){
    399                 if(!element->IsOnBase()) return;
    400                 basalelement=element->SpawnBasalElement();
    401         }
    402         else{
    403                 basalelement = element;
    404         }
     408
    405409        /*Fetch number of nodes and dof for this finite element*/
    406         int numnodes = basalelement->GetNumberOfNodes();
     410        int numnodes = element->GetNumberOfNodes();
    407411
    408412        /*Fetch dof list and allocate solution vector*/
    409         basalelement->GetDofList(&doflist,NoneApproximationEnum,GsetEnum);
     413        element->GetDofList(&doflist,NoneApproximationEnum,GsetEnum);
    410414        IssmDouble* newdamage = xNew<IssmDouble>(numnodes);
    411415
    412416        /*Get user-supplied max_damage: */
    413         basalelement->FindParam(&max_damage,DamageMaxDamageEnum);
     417        element->FindParam(&max_damage,DamageMaxDamageEnum);
    414418
    415419        /*Use the dof list to index into the solution vector: */
    … …  
    425429        /*Get all inputs and parameters*/
    426430        if(domaintype==Domain2DhorizontalEnum){
    427                 basalelement->AddInput(DamageDbarEnum,newdamage,element->GetElementType());
     431                element->AddInput(DamageDbarEnum,newdamage,element->GetElementType());
    428432        }
    429433        else{
    430                 basalelement->AddBasalInput(DamageDEnum,newdamage,element->GetElementType());
     434                element->AddInput(DamageDEnum,newdamage,element->GetElementType());
    431435        }
    432436
    … …  
    434438        xDelete<IssmDouble>(newdamage);
    435439        xDelete<int>(doflist);
    436         if(domaintype!=Domain2DhorizontalEnum){basalelement->DeleteMaterials(); delete basalelement;};
    437440}/*}}}*/
    438441void DamageEvolutionAnalysis::UpdateConstraints(FemModel* femmodel){/*{{{*/
    … …  
    446449        /*Intermediaries */
    447450        IssmDouble c1,c2,c3,healing,stress_threshold;
    448         IssmDouble s_xx,s_xy,s_yy,s1,s2,stmp;
     451        IssmDouble s_xx,s_xy,s_xz,s_yy,s_yz,s_zz,s1,s2,s3,stmp;
    449452        IssmDouble J2s,Chi,Psi,PosPsi,NegPsi;
    450         IssmDouble damage,tau_xx,tau_xy,tau_yy;
    451         int equivstress,domaintype,damagelaw;
     453        IssmDouble damage,tau_xx,tau_xy,tau_xz,tau_yy,tau_yz,tau_zz,stressMaxPrincipal;
     454        int equivstress,domaintype,damagelaw,dim;
    452455
    453456        /*Fetch number of vertices and allocate output*/
    … …  
    464467        element->FindParam(&damagelaw,DamageLawEnum);
    465468
    466         /*Compute stress tensor: */
     469        /*Get problem dimension*/
     470        switch(domaintype){
     471                case Domain2DhorizontalEnum: dim = 2; break;
     472                case Domain3DEnum:           dim = 3; break;
     473                default: _error_("not implemented");
     474        }
     475        /*Compute stress tensor and Stress Max Principal: */
    467476        element->ComputeDeviatoricStressTensor();
    468 
     477        if(dim==3){
     478                /*Only works in 3d because the pressure is defined*/
     479                element->StressMaxPrincipalCreateInput();
     480        }
    469481        /*retrieve what we need: */
    470482        Input* tau_xx_input  = element->GetInput(DeviatoricStressxxEnum);     _assert_(tau_xx_input);
    471483        Input* tau_xy_input  = element->GetInput(DeviatoricStressxyEnum);     _assert_(tau_xy_input);
    472484        Input* tau_yy_input  = element->GetInput(DeviatoricStressyyEnum);     _assert_(tau_yy_input);
     485        Input* tau_xz_input  = NULL;
     486        Input* tau_yz_input  = NULL;
     487        Input* tau_zz_input  = NULL;
     488        Input* stressMaxPrincipal_input = NULL;
     489        if(dim==3){
     490                tau_xz_input  = element->GetInput(DeviatoricStressxzEnum);     _assert_(tau_xz_input);
     491                tau_yz_input  = element->GetInput(DeviatoricStressyzEnum);     _assert_(tau_yz_input);
     492                tau_zz_input  = element->GetInput(DeviatoricStresszzEnum);     _assert_(tau_zz_input);
     493                stressMaxPrincipal_input = element->GetInput(StressMaxPrincipalEnum); _assert_(stressMaxPrincipal_input);
     494        }
    473495        Input* damage_input = NULL;
    474496        if(domaintype==Domain2DhorizontalEnum){
    … …  
    491513                tau_xy_input->GetInputValue(&tau_xy,gauss);
    492514                tau_yy_input->GetInputValue(&tau_yy,gauss);
    493        
     515                if(dim==3){
     516                        tau_xz_input->GetInputValue(&tau_xz,gauss);
     517                        tau_yz_input->GetInputValue(&tau_yz,gauss);
     518                        tau_zz_input->GetInputValue(&tau_zz,gauss);
     519                }
    494520                /*Calculate effective stress components*/
    495521                s_xx=tau_xx/(1.-damage);
    496522                s_xy=tau_xy/(1.-damage);
    497523                s_yy=tau_yy/(1.-damage);
    498 
     524                if(dim==3){
     525                        s_xz=tau_xz/(1.-damage);
     526                        s_yz=tau_yz/(1.-damage);
     527                        s_zz=tau_zz/(1.-damage);
     528                }
    499529                /*Calculate principal effective stresses*/
    500                 s1=(s_xx+s_yy)/2.+sqrt(pow((s_xx-s_yy)/2.,2)+pow(s_xy,2));
    501                 s2=(s_xx+s_yy)/2.-sqrt(pow((s_xx-s_yy)/2.,2)+pow(s_xy,2));
    502                 if(fabs(s2)>fabs(s1)){stmp=s2; s2=s1; s1=stmp;}
    503 
    504                 if(equivstress==0){ /* von Mises */
    505                         Chi=sqrt(s1*s1-s1*s2+s2*s2);
    506                 }
    507                 else if(equivstress==1){ /* max principal stress */
    508                         Chi=s1;
    509                 }
    510                 Psi=Chi-stress_threshold;
    511                 NegPsi=max(-Chi,0.); /* healing only for compressive stresses */
    512 
    513                 if(damagelaw==1){
    514                         PosPsi=max(Psi,0.);
    515                         f[i]= c1*(pow(PosPsi,c2) - healing*pow(NegPsi,c2))*pow((1./(1.-damage)),c3);
    516                 }
    517                 else if(damagelaw==2){
    518                         PosPsi=max(Psi,1.);
    519                         f[i]= c1*(pow(log10(PosPsi),c2) - healing*pow(NegPsi,c2))*pow((1./(1.-damage)),c3);
    520                 }
    521                 else _error_("damage law not supported");
    522         }
    523 
     530                if(dim==2){
     531                        s1=(s_xx+s_yy)/2.+sqrt(pow((s_xx-s_yy)/2.,2)+pow(s_xy,2));
     532                        s2=(s_xx+s_yy)/2.-sqrt(pow((s_xx-s_yy)/2.,2)+pow(s_xy,2));
     533                        if(fabs(s2)>fabs(s1)){stmp=s2; s2=s1; s1=stmp;}
     534
     535                        if(equivstress==0){ /* von Mises */
     536                                Chi=sqrt(s1*s1-s1*s2+s2*s2);
     537                        }
     538                        else if(equivstress==1){ /* max principal stress */
     539                                Chi=s1;
     540                        }
     541                        Psi=Chi-stress_threshold;
     542                        NegPsi=max(-Chi,0.); /* healing only for compressive stresses */
     543
     544                        if(damagelaw==1){
     545                                PosPsi=max(Psi,0.);
     546                                f[i]= c1*(pow(PosPsi,c2) - healing*pow(NegPsi,c2))*pow((1./(1.-damage)),c3);
     547                        }
     548                        else if(damagelaw==2){
     549                                PosPsi=max(Psi,1.);
     550                                f[i]= c1*(pow(log10(PosPsi),c2) - healing*pow(NegPsi,c2))*pow((1./(1.-damage)),c3);
     551                        }
     552                        else _error_("damage law not supported");
     553                }
     554                else{
     555                        if(equivstress==1){/* max principal stress */
     556                                stressMaxPrincipal_input->GetInputValue(&stressMaxPrincipal,gauss);
     557                                Chi=stressMaxPrincipal/(1.-damage);
     558                        }
     559                        else if(equivstress==0){/* von Mises */
     560                                Chi=sqrt(((s_xx-s_yy)*(s_xx-s_yy)+(s_yy-s_zz)*(s_yy-s_zz)+(s_zz-s_xx)*(s_zz-s_xx)+6.*(s_xy*s_xy+s_yz*s_yz+s_xz*s_xz))/2.);
     561                        }
     562                        Psi=Chi-stress_threshold;
     563                        NegPsi=max(-Chi,0.); /* healing only for compressive stresses */
     564                        if(damagelaw==1){
     565                                PosPsi=max(Psi,0.);
     566                                f[i]= c1*(pow(PosPsi,c2) - healing*pow(NegPsi,c2))*pow((1./(1.-damage)),c3);
     567                        }
     568                        else if(damagelaw==2){
     569                                PosPsi=max(Psi,1.);
     570                                f[i]= c1*(pow(log10(PosPsi),c2) - healing*pow(NegPsi,c2))*pow((1./(1.-damage)),c3);
     571                        }
     572                        else _error_("damage law not supported");
     573                }
     574        }
    524575        /*Add input*/
    525576        element->AddInput(DamageFEnum,f,element->GetElementType());

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

    r18403 r18470  
    14081408                /*Initialize maximum eigne value*/
    14091409                if(numroots>0){
    1410                         max = x[0];
     1410                        max = fabs(x[0]);
    14111411                }
    14121412                else{
    … …  
    14161416                /*Get max*/
    14171417                for(int i=1;i<numroots;i++){
    1418                         if(x[i]>max) max = x[i];
     1418                        if(fabs(x[i])>max) max = fabs(x[i]);
    14191419                }
    14201420

  • issm/trunk-jpl/src/c/classes/Node.cpp

    r18237 r18470  
    9494                                analysis_enum==BalancethicknessAnalysisEnum ||
    9595                                analysis_enum==HydrologyDCInefficientAnalysisEnum ||
    96                                 analysis_enum==DamageEvolutionAnalysisEnum ||
    9796                                analysis_enum==HydrologyDCEfficientAnalysisEnum ||
    9897                                analysis_enum==LevelsetAnalysisEnum ||

  • issm/trunk-jpl/src/c/cores/damage_core.cpp

    r18237 r18470  
    3434                femmodel->RequestedOutputsx(&femmodel->results,requested_outputs,numoutputs);
    3535        }
    36        
     36
    3737        /*Free resources:*/     
    3838        if(numoutputs){
Note: See TracChangeset for help on using the changeset viewer.