Changeset 16188


Ignore:
Timestamp:
09/19/13 17:33:08 (12 years ago)
Author:
Eric.Larour
Message:

CHG: finished putting together damage evolution framework. Need validation, verification.

Location:
issm/trunk-jpl
Files:
9 added
1 deleted
35 edited

Legend:

Unmodified
Added
Removed

  • issm/trunk-jpl/src/c/Makefile.am

    r16181 r16188  
    259259                                        ./modules/ModelProcessorx/CreateElementsVerticesAndMaterials.cpp\
    260260                                        ./modules/ModelProcessorx/CreateNodes.cpp\
     261                                        ./modules/ConstraintsStatex/PengridIsPresent.cpp\
    261262                                        ./modules/ParseToolkitsOptionsx/ParseToolkitsOptionsx.h\
    262263                                        ./modules/ParseToolkitsOptionsx/ParseToolkitsOptionsx.cpp\
    … …  
    409410                                           ./modules/ModelProcessorx/Melting/CreateLoadsMelting.cpp\
    410411                                           ./modules/ConstraintsStatex/ThermalConstraintsState.cpp\
    411                                            ./modules/ConstraintsStatex/ThermalIsPresent.cpp\
    412412                                           ./modules/ResetConstraintsx/ThermalConstraintsReset.cpp \
    413413                                           ./analyses/thermal_core.cpp\
    … …  
    551551                                        ./modules/ModelProcessorx/Damage/CreateNodesDamage.cpp \
    552552                                        ./modules/ModelProcessorx/Damage/CreateConstraintsDamage.cpp\
     553                                        ./modules/ConstraintsStatex/DamageConstraintsState.cpp\
     554                                        ./modules/ResetConstraintsx/DamageConstraintsReset.cpp \
    553555                                        ./modules/ModelProcessorx/Damage/CreateParametersDamage.cpp\
    554                                         ./modules/ModelProcessorx/Damage/CreateLoadsDamage.cpp
     556                                        ./modules/ModelProcessorx/Damage/CreateLoadsDamage.cpp\
     557                                        ./solutionsequences/solutionsequence_damage_nonlinear.cpp
    555558
    556559#}}}

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

    r16181 r16188  
    1 /*!\file: damage_core.cpp
    2  * \brief: core of the damage solution
     1/* 
     2 * \brief: damgage_core.cpp: core for the damage solution
    33 */
    44
    … …  
    1212void damage_core(FemModel* femmodel){
    1313
    14         /*parameters: */
    15         bool  save_results;
    16         bool  issmbgradients,ispdd,isdelta18o,isFS,isfreesurface;
    17         int   solution_type;
    18         int  *requested_outputs = NULL;
    19         int   numoutputs        = 0;
     14        /*intermediary*/
     15        bool   save_results;
     16        bool   dakota_analysis  = false;
     17        int    solution_type;
    2018
    21         /*activate configuration*/
    22         femmodel->SetCurrentConfiguration(DamageEvolutionAnalysisEnum);
    23          
    24         if(VerboseSolution())_printf_(" call damage core\n\n");
     19        if(VerboseSolution()) _printf0_("   computing damage\n");
     20       
     21        //first recover parameters common to all solutions
     22        femmodel->parameters->FindParam(&save_results,SaveResultsEnum);
     23        femmodel->parameters->FindParam(&dakota_analysis,QmuIsdakotaEnum);
     24        femmodel->parameters->FindParam(&solution_type,SolutionTypeEnum);
    2525
    26         /*recover parameters: */
    27         femmodel->parameters->FindParam(&save_results,SaveResultsEnum);
    28         femmodel->parameters->FindParam(&issmbgradients,SurfaceforcingsIssmbgradientsEnum);
    29         femmodel->parameters->FindParam(&ispdd,SurfaceforcingsIspddEnum);
    30         femmodel->parameters->FindParam(&isdelta18o,SurfaceforcingsIsdelta18oEnum);
    31         femmodel->parameters->FindParam(&isFS,FlowequationIsFSEnum);
    32         femmodel->parameters->FindParam(&isfreesurface,MasstransportIsfreesurfaceEnum);
    33         femmodel->parameters->FindParam(&solution_type,SolutionTypeEnum);
    34         femmodel->parameters->FindParam(&numoutputs,MasstransportNumRequestedOutputsEnum);
    35         if(numoutputs) femmodel->parameters->FindParam(&requested_outputs,&numoutputs,MasstransportRequestedOutputsEnum);
    36 
    37         if(issmbgradients){
    38           if(VerboseSolution())_printf_("       call smb gradients module\n\n");
    39           SmbGradientsx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters);
    40         }
    41         if(ispdd){
    42                 if(isdelta18o){
    43                         if(VerboseSolution()) _printf0_("   call Delta18oParametrization module\n");
    44                         Delta18oParameterizationx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters);
    45                 }
    46                 if(VerboseSolution()) _printf0_("   call positive degree day module\n");
    47                 PositiveDegreeDayx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters);
     26        if(dakota_analysis && solution_type!=TransientSolutionEnum){
     27                femmodel->SetCurrentConfiguration(DamageEvolutionAnalysisEnum);
     28                ResetConstraintsx(femmodel);
    4829        }
    4930
    50         if(isFS && isfreesurface){
    51                 if(VerboseSolution()) _printf0_("   call free surface computational core\n");
    52                 femmodel->SetCurrentConfiguration(FreeSurfaceBaseAnalysisEnum);
    53                 solutionsequence_linear(femmodel);
    54                 femmodel->SetCurrentConfiguration(FreeSurfaceTopAnalysisEnum);
    55                 solutionsequence_linear(femmodel);
    56         }
    57         else{
    58                 if(VerboseSolution()) _printf0_("   call computational core\n");
    59                 solutionsequence_linear(femmodel);
    60         }
     31        femmodel->SetCurrentConfiguration(DamageEvolutionAnalysisEnum);
     32        solutionsequence_damage_nonlinear(femmodel);
    6133
    6234        if(save_results){
    6335                if(VerboseSolution()) _printf0_("   saving results\n");
    64                 InputToResultx(femmodel,ThicknessEnum);
    65                 InputToResultx(femmodel,BedEnum);
    66                 InputToResultx(femmodel,SurfaceEnum);
    67                 femmodel->RequestedOutputsx(requested_outputs,numoutputs);
     36                InputToResultx(femmodel,DamageDEnum);
    6837        }
    69 
    70         if(solution_type==MasstransportSolutionEnum)femmodel->RequestedDependentsx();
    71 
    72         /*Free ressources:*/
    73         xDelete<int>(requested_outputs);
    7438}

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

    r16186 r16188  
    5454                virtual void   GetInputListOnVertices(IssmDouble* pvalue,int enumtype,IssmDouble defaultvalue)=0;
    5555                virtual void   GetInputValue(IssmDouble* pvalue,Node* node,int enumtype)=0;
     56                virtual void   GetMaterialInputValue(IssmDouble* pvalue,Node* node,int enumtype)=0;
    5657
    5758                virtual IssmDouble SurfaceArea(void)=0;

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

    r16186 r16188  
    14221422
    14231423        Input* input=inputs->GetInput(enumtype);
     1424        if(!input) _error_("No input of type " << EnumToStringx(enumtype) << " found in tria");
     1425
     1426        GaussPenta* gauss=new GaussPenta();
     1427        gauss->GaussVertex(this->GetNodeIndex(node));
     1428
     1429        input->GetInputValue(pvalue,gauss);
     1430        delete gauss;
     1431}
     1432/*}}}*/
     1433/*FUNCTION Penta::GetMaterialInputValue(IssmDouble* pvalue,Node* node,int enumtype) {{{*/
     1434void Penta::GetMaterialInputValue(IssmDouble* pvalue,Node* node,int enumtype){
     1435
     1436        Input* input=this->material->inputs->GetInput(enumtype);
    14241437        if(!input) _error_("No input of type " << EnumToStringx(enumtype) << " found in tria");
    14251438

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

    r16186 r16188  
    200200                void           GetInputListOnNodes(IssmDouble* pvalue,int enumtype,IssmDouble defaultvalue);
    201201                void           GetInputValue(IssmDouble* pvalue,Node* node,int enumtype);
     202                void           GetMaterialInputValue(IssmDouble* pvalue,Node* node,int enumtype);
    202203                void             GetPhi(IssmDouble* phi, IssmDouble*  epsilon, IssmDouble viscosity);
    203204                void           GetQuadNormal(IssmDouble* normal,IssmDouble xyz_list[4][3]);

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

    r16186 r16188  
    176176void  Tria::CreateKMatrix(Matrix<IssmDouble>* Kff, Matrix<IssmDouble>* Kfs){
    177177
     178
    178179        /*Skip if water element*/
    179180        if(NoIceInElement()) return;
    … …  
    223224                        return CreateKMatrixStressbalanceSIA();
    224225                        break;
    225                  #endif
     226                #endif
     227                #ifdef _HAVE_DAMAGE_
     228                case DamageEvolutionAnalysisEnum:
     229                        return CreateKMatrixDamageEvolution();
     230                        break;
     231                #endif
    226232                case BedSlopeXAnalysisEnum: case SurfaceSlopeXAnalysisEnum: case BedSlopeYAnalysisEnum: case SurfaceSlopeYAnalysisEnum:
    227233                        return CreateMassMatrix();
    … …  
    363369        /*Just branch to the correct load generator, according to the type of analysis we are carrying out: */
    364370        switch(analysis_type){
    365 #ifdef _HAVE_STRESSBALANCE_
     371                #ifdef _HAVE_STRESSBALANCE_
    366372                case StressbalanceAnalysisEnum:
    367373                        return CreatePVectorStressbalanceSSA();
    … …  
    370376                        return CreatePVectorStressbalanceSIA();
    371377                        break;
    372 #endif
     378                #endif
     379                #ifdef _HAVE_DAMAGE_
     380                case DamageEvolutionAnalysisEnum:
     381                        return CreatePVectorDamageEvolution();
     382                        break;
     383                #endif
    373384                case BedSlopeXAnalysisEnum: case SurfaceSlopeXAnalysisEnum: case BedSlopeYAnalysisEnum: case SurfaceSlopeYAnalysisEnum:
    374385                        return CreatePVectorSlope();
    … …  
    377388                        return CreatePVectorMasstransport();
    378389                        break;
    379 #ifdef _HAVE_HYDROLOGY_
     390                #ifdef _HAVE_HYDROLOGY_
    380391                case HydrologyShreveAnalysisEnum:
    381392                        return CreatePVectorHydrologyShreve();
    … …  
    387398                        return CreatePVectorHydrologyDCEfficient();
    388399                        break;
    389 #endif
    390 #ifdef _HAVE_BALANCED_
     400                #endif
     401                #ifdef _HAVE_BALANCED_
    391402                case BalancethicknessAnalysisEnum:
    392403                        return CreatePVectorBalancethickness();
    … …  
    401412                        return CreatePVectorSmoothedSlopeY();
    402413                        break;
    403 #endif
    404 #ifdef _HAVE_CONTROL_
     414                #endif
     415                #ifdef _HAVE_CONTROL_
    405416                case AdjointBalancethicknessAnalysisEnum:
    406417                        return CreatePVectorAdjointBalancethickness();
    … …  
    409420                        return CreatePVectorAdjointHoriz();
    410421                        break;
    411 #endif
     422                #endif
    412423                default:
    413424                        _error_("analysis " << analysis_type << " (" << EnumToStringx(analysis_type) << ") not supported yet");
    … …  
    12391250
    12401251        Input* input=inputs->GetInput(enumtype);
     1252        if(!input) _error_("No input of type " << EnumToStringx(enumtype) << " found in tria");
     1253
     1254        GaussTria* gauss=new GaussTria();
     1255        gauss->GaussVertex(this->GetNodeIndex(node));
     1256
     1257        input->GetInputValue(pvalue,gauss);
     1258        delete gauss;
     1259}
     1260/*}}}*/
     1261/*FUNCTION Tria::GetMaterialInputValue(IssmDouble* pvalue,Node* node,int enumtype) {{{*/
     1262void Tria::GetMaterialInputValue(IssmDouble* pvalue,Node* node,int enumtype){
     1263
     1264        Input* input=this->material->inputs->GetInput(enumtype);
    12411265        if(!input) _error_("No input of type " << EnumToStringx(enumtype) << " found in tria");
    12421266
    … …  
    16231647                case HydrologyDCEfficientAnalysisEnum:
    16241648                        InputUpdateFromSolutionOneDof(solution,EplHeadEnum);
     1649                        break;
     1650                #endif
     1651                #ifdef _HAVE_DAMAGE_
     1652                case DamageEvolutionAnalysisEnum:
     1653                        InputUpdateFromSolutionDamageEvolution(solution);
    16251654                        break;
    16261655                #endif
    … …  
    68256854#endif
    68266855
     6856#ifdef _HAVE_DAMAGE_
     6857/*FUNCTION Tria::CreateKMatrixDamageEvolution {{{*/
     6858ElementMatrix* Tria::CreateKMatrixDamageEvolution(void){
     6859
     6860        switch(GetElementType()){
     6861                case P1Enum: case P2Enum:
     6862                        return CreateKMatrixDamageEvolution_CG();
     6863                case P1DGEnum:
     6864                        _error_("DG not implemented yet!");break;
     6865                default:
     6866                        _error_("Element type " << EnumToStringx(GetElementType()) << " not supported yet");
     6867        }
     6868}
     6869/*}}}*/
     6870/*FUNCTION Tria::CreateKMatrixDamageEvolution_CG {{{*/
     6871ElementMatrix* Tria::CreateKMatrixDamageEvolution_CG(void){
     6872
     6873        /*Intermediaries */
     6874        int        stabilization;
     6875        int        dim;
     6876        IssmDouble Jdet,D_scalar,dt,h;
     6877        IssmDouble vel,vx,vy,dvxdx,dvydy;
     6878        IssmDouble dvx[2],dvy[2];
     6879        IssmDouble xyz_list[NUMVERTICES][3];
     6880
     6881        /*Fetch number of nodes for this finite element*/
     6882        int numnodes = this->NumberofNodes();
     6883
     6884        /*Initialize Element matrix and vectors*/
     6885        ElementMatrix* Ke     = new ElementMatrix(nodes,numnodes,this->parameters,NoneApproximationEnum);
     6886        IssmDouble*    basis  = xNew<IssmDouble>(numnodes);
     6887        IssmDouble*    B      = xNew<IssmDouble>(2*numnodes);
     6888        IssmDouble*    Bprime = xNew<IssmDouble>(2*numnodes);
     6889        IssmDouble     D[2][2];
     6890
     6891        /*Retrieve all inputs and parameters*/
     6892        GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
     6893        this->parameters->FindParam(&dt,TimesteppingTimeStepEnum);
     6894        this->parameters->FindParam(&dim,MeshDimensionEnum);
     6895        this->parameters->FindParam(&stabilization,DamageStabilizationEnum);
     6896        Input* vxaverage_input=NULL;
     6897        Input* vyaverage_input=NULL;
     6898        if(dim==2){
     6899                vxaverage_input=inputs->GetInput(VxEnum); _assert_(vxaverage_input);
     6900                vyaverage_input=inputs->GetInput(VyEnum); _assert_(vyaverage_input);
     6901        }
     6902        else{
     6903                vxaverage_input=inputs->GetInput(VxAverageEnum); _assert_(vxaverage_input);
     6904                vyaverage_input=inputs->GetInput(VyAverageEnum); _assert_(vyaverage_input);
     6905        }
     6906        h=sqrt(2*this->GetArea());
     6907
     6908        /* Start  looping on the number of gaussian points: */
     6909        GaussTria *gauss=new GaussTria(2);
     6910        for(int ig=gauss->begin();ig<gauss->end();ig++){
     6911
     6912                gauss->GaussPoint(ig);
     6913
     6914                GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss);
     6915                GetNodalFunctions(basis,gauss);
     6916
     6917                vxaverage_input->GetInputValue(&vx,gauss);
     6918                vyaverage_input->GetInputValue(&vy,gauss);
     6919                vxaverage_input->GetInputDerivativeValue(&dvx[0],&xyz_list[0][0],gauss);
     6920                vyaverage_input->GetInputDerivativeValue(&dvy[0],&xyz_list[0][0],gauss);
     6921
     6922                D_scalar=gauss->weight*Jdet;
     6923
     6924                TripleMultiply(basis,1,numnodes,1,
     6925                                        &D_scalar,1,1,0,
     6926                                        basis,1,numnodes,0,
     6927                                        &Ke->values[0],1);
     6928                GetBMasstransport(B,&xyz_list[0][0],gauss);
     6929                GetBprimeMasstransport(Bprime,&xyz_list[0][0],gauss);
     6930
     6931                dvxdx=dvx[0];
     6932                dvydy=dvy[1];
     6933                D_scalar=dt*gauss->weight*Jdet;
     6934
     6935                D[0][0]=D_scalar*dvxdx;
     6936                D[0][1]=0.;
     6937                D[1][0]=0.;
     6938                D[1][1]=D_scalar*dvydy;
     6939                TripleMultiply(B,2,numnodes,1,
     6940                                        &D[0][0],2,2,0,
     6941                                        B,2,numnodes,0,
     6942                                        &Ke->values[0],1);
     6943
     6944                D[0][0]=D_scalar*vx;
     6945                D[1][1]=D_scalar*vy;
     6946                TripleMultiply(B,2,numnodes,1,
     6947                                        &D[0][0],2,2,0,
     6948                                        Bprime,2,numnodes,0,
     6949                                        &Ke->values[0],1);
     6950
     6951                if(stabilization==2){
     6952                        /*Streamline upwinding*/
     6953                        vel=sqrt(vx*vx+vy*vy)+1.e-8;
     6954                        D[0][0]=h/(2*vel)*vx*vx;
     6955                        D[1][0]=h/(2*vel)*vy*vx;
     6956                        D[0][1]=h/(2*vel)*vx*vy;
     6957                        D[1][1]=h/(2*vel)*vy*vy;
     6958                }
     6959                else if(stabilization==1){
     6960                        /*SSA*/
     6961                        vxaverage_input->GetInputAverage(&vx);
     6962                        vyaverage_input->GetInputAverage(&vy);
     6963                        D[0][0]=h/2.0*fabs(vx);
     6964                        D[0][1]=0.;
     6965                        D[1][0]=0.;
     6966                        D[1][1]=h/2.0*fabs(vy);
     6967                }
     6968                if(stabilization==1 || stabilization==2){
     6969                        D[0][0]=D_scalar*D[0][0];
     6970                        D[1][0]=D_scalar*D[1][0];
     6971                        D[0][1]=D_scalar*D[0][1];
     6972                        D[1][1]=D_scalar*D[1][1];
     6973                        TripleMultiply(Bprime,2,numnodes,1,
     6974                                                &D[0][0],2,2,0,
     6975                                                Bprime,2,numnodes,0,
     6976                                                &Ke->values[0],1);
     6977                }
     6978        }
     6979
     6980        /*Clean up and return*/
     6981        xDelete<IssmDouble>(basis);
     6982        xDelete<IssmDouble>(B);
     6983        xDelete<IssmDouble>(Bprime);
     6984        delete gauss;
     6985        return Ke;
     6986}
     6987/*}}}*/
     6988/*FUNCTION Tria::CreatePVectorDamageEvolution{{{*/
     6989ElementVector* Tria::CreatePVectorDamageEvolution(void){
     6990
     6991        switch(GetElementType()){
     6992                case P1Enum: case P2Enum:
     6993                        return CreatePVectorDamageEvolution_CG();
     6994                case P1DGEnum:
     6995                        _error_("DG not implemented yet");
     6996                default:
     6997                        _error_("Element type " << EnumToStringx(GetElementType()) << " not supported yet");
     6998        }
     6999}
     7000/*}}}*/
     7001/*FUNCTION Tria::CreatePVectorDamageEvolution_CG {{{*/
     7002ElementVector* Tria::CreatePVectorDamageEvolution_CG(void){
     7003
     7004        /*Intermediaries */
     7005        IssmDouble Jdet,dt;
     7006        IssmDouble f,damage;
     7007        IssmDouble xyz_list[NUMVERTICES][3];
     7008        Input* damage_input  = NULL;
     7009        Input* f_input  = NULL;
     7010
     7011        /*Fetch number of nodes and dof for this finite element*/
     7012        int numnodes = this->NumberofNodes();
     7013
     7014        /*Initialize Element vector and other vectors*/
     7015        ElementVector* pe    = new ElementVector(nodes,numnodes,this->parameters);
     7016        IssmDouble*    basis = xNew<IssmDouble>(numnodes);
     7017
     7018        /*Retrieve all inputs and parameters*/
     7019        GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
     7020        this->parameters->FindParam(&dt,TimesteppingTimeStepEnum);
     7021        damage_input  = this->material->inputs->GetInput(DamageDbarEnum);     _assert_(damage_input);
     7022       
     7023        /*retrieve damage evolution forcing function: */
     7024        f_input=this->DamageEvolutionF();
     7025
     7026        /*Initialize forcing function f to 0, do not forget!:*/
     7027        /* Start  looping on the number of gaussian points: */
     7028        GaussTria* gauss=new GaussTria(2);
     7029        for(int ig=gauss->begin();ig<gauss->end();ig++){
     7030
     7031                gauss->GaussPoint(ig);
     7032
     7033                GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss);
     7034                GetNodalFunctions(basis,gauss);
     7035
     7036                f_input->GetInputValue(&f,gauss);
     7037                damage_input->GetInputValue(&damage,gauss);
     7038
     7039                for(int i=0;i<numnodes;i++) pe->values[i]+=Jdet*gauss->weight*(damage+dt*f)*basis[i];
     7040        }
     7041
     7042        /*Clean up and return*/
     7043        xDelete<IssmDouble>(basis);
     7044        delete gauss;
     7045        return pe;
     7046}
     7047/*}}}*/
     7048/*FUNCTION Tria::DamageEvolutionF{{{*/
     7049Input* Tria::DamageEvolutionF(void){
     7050
     7051        /*Intermediaries */
     7052        IssmDouble c1,c2,c3,healing,stress_threshold;
     7053        IssmDouble s_xx,s_xy,s_yy;
     7054        IssmDouble J2s;
     7055        IssmDouble Xis;
     7056        IssmDouble Psi;
     7057        IssmDouble PosPsi;
     7058        IssmDouble NegPsi;
     7059        Input* damage_input=NULL;
     7060        Input* z_input=NULL;
     7061        Input* sigma_xx_input  = NULL;
     7062        Input* sigma_xy_input  = NULL;
     7063        Input* sigma_yy_input  = NULL;
     7064        GaussTria* gauss=NULL;
     7065        IssmDouble damage,sigma_xx,sigma_xy,sigma_yy;
     7066       
     7067        /*output: */
     7068        IssmDouble f[NUMVERTICES];
     7069        Input* f_input=NULL;
     7070
     7071        /*retrieve parameters:*/
     7072        this->parameters->FindParam(&c1,DamageC1Enum);
     7073        this->parameters->FindParam(&c2,DamageC2Enum);
     7074        this->parameters->FindParam(&c3,DamageC3Enum);
     7075        this->parameters->FindParam(&healing,DamageHealingEnum);
     7076        this->parameters->FindParam(&stress_threshold,DamageStressThresholdEnum);
     7077
     7078        /*Compute stress tensor: */
     7079        this->ComputeStressTensor();
     7080
     7081        /*retrieve what we need: */
     7082        sigma_xx_input  = inputs->GetInput(StressTensorxxEnum);  _assert_(sigma_xx_input);
     7083        sigma_xy_input  = inputs->GetInput(StressTensorxyEnum);  _assert_(sigma_xy_input);
     7084        sigma_yy_input  = inputs->GetInput(StressTensoryyEnum);  _assert_(sigma_yy_input);
     7085        damage_input  = this->material->inputs->GetInput(DamageDbarEnum);        _assert_(damage_input);
     7086
     7087        /*Damage evolution z mapping: */
     7088        gauss=new GaussTria();
     7089        J2s=0;
     7090        for (int iv=0;iv<NUMVERTICES;iv++){
     7091                gauss->GaussVertex(iv);
     7092               
     7093                damage_input->GetInputValue(&damage,gauss);
     7094                sigma_xx_input->GetInputValue(&sigma_xx,gauss);
     7095                sigma_xy_input->GetInputValue(&sigma_xy,gauss);
     7096                sigma_yy_input->GetInputValue(&sigma_yy,gauss);
     7097
     7098                s_xx=sigma_xx/(1-damage);
     7099                s_xy=sigma_xy/(1-damage);
     7100                s_yy=sigma_yy/(1-damage);
     7101
     7102                J2s=1.0/sqrt(2.0)*sqrt(pow(s_xx,2)+pow(s_yy,2)+2*pow(s_xy,2));
     7103               
     7104                Xis=sqrt(3.0)*J2s;
     7105
     7106                Psi=Xis-stress_threshold;
     7107
     7108                PosPsi=max(Psi,0.0);
     7109                NegPsi=max(-Psi,0.0);
     7110
     7111                f[iv]= c1* ( pow(PosPsi,c2)   -  healing * pow(NegPsi,c2) )  *   pow((1 - damage),-c3);
     7112
     7113        }
     7114        f_input=new TriaInput(NoneEnum,&f[0],P1Enum);
     7115       
     7116        /*Clean up and return*/
     7117        delete gauss;
     7118        return f_input;
     7119}
     7120/*}}}*/
     7121/*FUNCTION Tria::InputUpdateFromSolutionDamageEvolution {{{*/
     7122void  Tria::InputUpdateFromSolutionDamageEvolution(IssmDouble* solution){
     7123
     7124        const int    numdof=NDOF1*NUMVERTICES;
     7125
     7126        int         i;
     7127        IssmDouble  values[numdof];
     7128        int        *doflist = NULL;
     7129
     7130        /*Get dof list: */
     7131        GetDofList(&doflist,NoneApproximationEnum,GsetEnum);
     7132
     7133        /*Use the dof list to index into the solution vector: */
     7134        for(i=0;i<numdof;i++){
     7135                values[i]=solution[doflist[i]];
     7136                /*Check solution*/
     7137                if(xIsNan<IssmDouble>(values[i])) _error_("NaN found in solution vector");
     7138        }
     7139
     7140        /*Get all inputs and parameters*/
     7141        this->material->inputs->AddInput(new TriaInput(DamageDbarEnum,values,P1Enum));
     7142
     7143        /*Free ressources:*/
     7144        xDelete<int>(doflist);
     7145}
     7146/*}}}*/
     7147#endif
     7148
    68277149#ifdef _HAVE_DAKOTA_
    68287150/*FUNCTION Tria::InputUpdateFromVectorDakota(IssmDouble* vector, int name, int type);{{{*/

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

    r16186 r16188  
    218218                void           GetInputListOnNodes(IssmDouble* pvalue,int enumtype,IssmDouble defaultvalue);
    219219                void           GetInputValue(IssmDouble* pvalue,Node* node,int enumtype);
     220                void           GetMaterialInputValue(IssmDouble* pvalue,Node* node,int enumtype);
    220221                void           GetStrainRate2d(IssmDouble* epsilon,IssmDouble* xyz_list, GaussTria* gauss, Input* vx_input, Input* vy_input);
    221222                void             InputUpdateFromSolutionOneDof(IssmDouble* solution,int enum_type);
    … …  
    276277                bool    AnyActive(void);
    277278                #endif
    278                 #ifdef _HAVE_BALANCED_
    279                 #endif
     279
     280                #ifdef _HAVE_DAMAGE_
     281                ElementMatrix* CreateKMatrixDamageEvolution(void);
     282                ElementMatrix* CreateKMatrixDamageEvolution_CG(void);
     283                ElementVector* CreatePVectorDamageEvolution(void);
     284                ElementVector* CreatePVectorDamageEvolution_CG(void);
     285                Input* DamageEvolutionF(void);
     286                void           InputUpdateFromSolutionDamageEvolution(IssmDouble* solution);
     287                #endif
     288
    280289
    281290                /*}}}*/

  • issm/trunk-jpl/src/c/classes/Inputs/BoolInput.h

    r15737 r16188  
    7272                IssmDouble Min(void){_error_("Min not implemented for booleans");};
    7373                IssmDouble MinAbs(void){_error_("Min not implemented for booleans");};
     74                void Set(IssmDouble setvalue){_error_("Set not implemented yet");};
    7475                void Scale(IssmDouble scale_factor);
    7576                void ArtificialNoise(IssmDouble min,IssmDouble max){_error_("not implemented yet");};

  • issm/trunk-jpl/src/c/classes/Inputs/ControlInput.h

    r15737 r16188  
    7272                void SquareMin(IssmDouble* psquaremin,Parameters* parameters){_error_("not implemented yet");};
    7373                void ConstrainMin(IssmDouble minimum){_error_("not implemented yet");};
     74                void Set(IssmDouble setvalue){_error_("Set not implemented yet");};
    7475                void Scale(IssmDouble scale_factor){_error_("not implemented yet");};
    7576                void ArtificialNoise(IssmDouble min,IssmDouble max){_error_("not implemented yet");};

  • issm/trunk-jpl/src/c/classes/Inputs/DatasetInput.h

    r15737 r16188  
    6666                void SquareMin(IssmDouble* psquaremin,Parameters* parameters){_error_("not implemented yet");};
    6767                void ConstrainMin(IssmDouble minimum){_error_("not implemented yet");};
     68                void Set(IssmDouble setvalue){_error_("Set not implemented yet");};
    6869                void Scale(IssmDouble scale_factor){_error_("not implemented yet");};
    6970                void ArtificialNoise(IssmDouble min,IssmDouble max){_error_("not implemented yet");};

  • issm/trunk-jpl/src/c/classes/Inputs/DoubleInput.h

    r15737 r16188  
    6666                void SquareMin(IssmDouble* psquaremin,Parameters* parameters);
    6767                void ConstrainMin(IssmDouble minimum);
     68                void Set(IssmDouble setvalue){_error_("Set not implemented yet");};
    6869                void Scale(IssmDouble scale_factor);
    6970                void ArtificialNoise(IssmDouble min,IssmDouble max){_error_("not implemented yet");};

  • issm/trunk-jpl/src/c/classes/Inputs/Input.h

    r15737 r16188  
    5656                virtual IssmDouble Max(void)=0;
    5757                virtual IssmDouble Min(void)=0;
     58                virtual void Set(IssmDouble setvalue)=0;
    5859                virtual void   Scale(IssmDouble scale_factor)=0;
    5960                virtual void   ArtificialNoise(IssmDouble min,IssmDouble max)=0;

  • issm/trunk-jpl/src/c/classes/Inputs/IntInput.h

    r15737 r16188  
    6868                void SquareMin(IssmDouble* psquaremin,Parameters* parameters);
    6969                void ConstrainMin(IssmDouble minimum){_error_("not implemented yet");};
     70                void Set(IssmDouble setvalue){_error_("Set not implemented yet");};
    7071                void Scale(IssmDouble scale_factor);
    7172                void ArtificialNoise(IssmDouble min,IssmDouble max){_error_("not implemented yet");};

  • issm/trunk-jpl/src/c/classes/Inputs/PentaInput.h

    r15737 r16188  
    6868                void SquareMin(IssmDouble* psquaremin,Parameters* parameters);
    6969                void ConstrainMin(IssmDouble minimum);
     70                void Set(IssmDouble setvalue){_error_("Set not implemented yet");};
    7071                void Scale(IssmDouble scale_factor);
    7172                void ArtificialNoise(IssmDouble min,IssmDouble max){_error_("not implemented yet");};

  • issm/trunk-jpl/src/c/classes/Inputs/TransientInput.h

    r15737 r16188  
    7171                void SquareMin(IssmDouble* psquaremin,Parameters* parameters);
    7272                void ConstrainMin(IssmDouble minimum){_error_("not implemented yet");};
     73                void Set(IssmDouble setvalue){_error_("Set not implemented yet");};
    7374                void Scale(IssmDouble scale_factor){_error_("not implemented yet");};
    7475                void ArtificialNoise(IssmDouble min,IssmDouble max){_error_("not implemented yet");};

  • issm/trunk-jpl/src/c/classes/Inputs/TriaInput.cpp

    r15737 r16188  
    321321}
    322322/*}}}*/
     323/*FUNCTION TriaInput::Set{{{*/
     324void TriaInput::Set(IssmDouble setvalue){
     325
     326        const int numnodes=this->NumberofNodes();
     327        for(int i=0;i<numnodes;i++)values[i]=setvalue;
     328}
     329/*}}}*/
    323330/*FUNCTION TriaInput::ArtificialNoise{{{*/
    324331void TriaInput::ArtificialNoise(IssmDouble min,IssmDouble max){
    … …  
    435442}
    436443/*}}}*/
     444/*FUNCTION TriaInput::PointwiseDivide{{{*/
     445Input* TriaInput::PointwiseDivide(Input* inputB){
     446
     447        /*Ouput*/
     448        TriaInput* outinput=NULL;
     449
     450        /*Intermediaries*/
     451        TriaInput *xinputB  = NULL;
     452        const int   numnodes = this->NumberofNodes();
     453
     454        /*Check that inputB is of the same type*/
     455        if(inputB->ObjectEnum()!=TriaInputEnum)     _error_("Operation not permitted because inputB is of type " << EnumToStringx(inputB->ObjectEnum()));
     456        xinputB=(TriaInput*)inputB;
     457        if(xinputB->element_type!=this->element_type) _error_("Operation not permitted because inputB is of type " << EnumToStringx(xinputB->element_type));
     458
     459        /*Allocate intermediary*/
     460        IssmDouble* AdotBvalues=xNew<IssmDouble>(numnodes);
     461
     462        /*Create point wise division*/
     463        for(int i=0;i<numnodes;i++){
     464                _assert_(xinputB->values[i]!=0);
     465                AdotBvalues[i]=this->values[i]/xinputB->values[i];
     466        }
     467
     468        /*Create new Tria vertex input (copy of current input)*/
     469        outinput=new TriaInput(this->enum_type,AdotBvalues,this->element_type);
     470
     471        /*Return output pointer*/
     472        xDelete<IssmDouble>(AdotBvalues);
     473        return outinput;
     474
     475}
     476/*}}}*/
    437477/*FUNCTION TriaInput::Configure{{{*/
    438478void TriaInput::Configure(Parameters* parameters){

  • issm/trunk-jpl/src/c/classes/Inputs/TriaInput.h

    r15737 r16188  
    3535                int    InstanceEnum();
    3636                Input* SpawnTriaInput(int location);
    37                 Input* PointwiseDivide(Input* inputB){_error_("not implemented yet");};
     37                Input* PointwiseDivide(Input* inputB);
    3838                Input* PointwiseMin(Input* inputB);
    3939                Input* PointwiseMax(Input* inputB);
    … …  
    6868                void SquareMin(IssmDouble* psquaremin,Parameters* parameters);
    6969                void ConstrainMin(IssmDouble minimum);
     70                void Set(IssmDouble setvalue);
    7071                void Scale(IssmDouble scale_factor);
    7172                void ArtificialNoise(IssmDouble min,IssmDouble max);

  • issm/trunk-jpl/src/c/classes/Loads/Pengrid.cpp

    r16166 r16188  
    241241                case HydrologyDCInefficientAnalysisEnum:
    242242                        Ke=PenaltyCreateKMatrixHydrologyDCInefficient(kmax);
     243                        break;
     244                #endif
     245                #ifdef _HAVE_DAMAGE_
     246                case DamageEvolutionAnalysisEnum:
     247                        Ke=PenaltyCreateKMatrixDamageEvolution(kmax);
    243248                        break;
    244249                #endif
    … …  
    278283                        break;
    279284                #endif
     285                #ifdef _HAVE_DAMAGE_
     286                case DamageEvolutionAnalysisEnum:
     287                        pe=PenaltyCreatePVectorDamageEvolution(kmax);
     288                        break;
     289                #endif
     290
    280291                default:
    281292                        _error_("analysis " << analysis_type << " (" << EnumToStringx(analysis_type) << ") not supported yet");
    … …  
    430441        }
    431442        #endif
     443        #ifdef _HAVE_DAMAGE_
     444        else if (analysis_type==DamageEvolutionAnalysisEnum){
     445                ConstraintActivateDamageEvolution(punstable);
     446                return;
     447        }
     448        #endif
     449
    432450        else{
    433451                _error_("analysis: " << EnumToStringx(analysis_type) << " not supported yet");
    … …  
    695713/*}}}*/
    696714#endif
     715#ifdef _HAVE_DAMAGE_
     716/*FUNCTION Pengrid::ConstraintActivateDamageEvolution {{{*/
     717void  Pengrid::ConstraintActivateDamageEvolution(int* punstable){
     718
     719        //   The penalty is stable if it doesn't change during to successive iterations.   
     720        IssmDouble max_damage;
     721        IssmDouble damage;
     722        int        new_active;
     723        int        unstable=0;
     724        int        penalty_lock;
     725
     726        /*check that pengrid is not a clone (penalty to be added only once)*/
     727        if (node->IsClone()){
     728                unstable=0;
     729                *punstable=unstable;
     730                return;
     731        }
     732
     733        //First recover damage  using the element: */
     734        element->GetMaterialInputValue(&damage,node,DamageDbarEnum);
     735
     736        //Recover our data:
     737        parameters->FindParam(&penalty_lock,DamagePenaltyLockEnum);
     738        parameters->FindParam(&max_damage,DamageMaxDamageEnum);
     739       
     740        //Figure out if damage is over max_damage, in which case, this penalty needs to be activated.
     741
     742        if (damage>max_damage){
     743                new_active=1;
     744        }
     745        else{
     746                new_active=0;
     747        }
     748
     749        //Figure out stability of this penalty
     750        if (active==new_active){
     751                unstable=0;
     752        }
     753        else{
     754                unstable=1;
     755                if(penalty_lock)zigzag_counter++;
     756        }
     757
     758        /*If penalty keeps zigzagging more than penalty_lock times: */
     759        if(penalty_lock){
     760                if(zigzag_counter>penalty_lock){
     761                        unstable=0;
     762                        active=1;
     763                }
     764        }
     765
     766        //Set penalty flag
     767        active=new_active;
     768
     769        //*Assign output pointers:*/
     770        *punstable=unstable;
     771}
     772/*}}}*/
     773/*FUNCTION Pengrid::PenaltyCreateKMatrixDamageEvolution {{{*/
     774ElementMatrix* Pengrid::PenaltyCreateKMatrixDamageEvolution(IssmDouble kmax){
     775
     776        IssmDouble    penalty_factor;
     777
     778        /*Initialize Element matrix and return if necessary*/
     779        if(!this->active) return NULL;
     780        ElementMatrix* Ke=new ElementMatrix(&node,NUMVERTICES,this->parameters);
     781
     782        /*recover parameters: */
     783        parameters->FindParam(&penalty_factor,DamagePenaltyFactorEnum);
     784
     785        Ke->values[0]=kmax*pow(10.,penalty_factor);
     786
     787        /*Clean up and return*/
     788        return Ke;
     789}
     790/*}}}*/
     791/*FUNCTION Pengrid::PenaltyCreatePVectorDamageEvolution {{{*/
     792ElementVector* Pengrid::PenaltyCreatePVectorDamageEvolution(IssmDouble kmax){
     793
     794        IssmDouble penalty_factor;
     795        IssmDouble max_damage;
     796
     797        /*Initialize Element matrix and return if necessary*/
     798        if(!this->active) return NULL;
     799        ElementVector* pe=new ElementVector(&node,1,this->parameters);
     800
     801        //Recover our data:
     802        parameters->FindParam(&penalty_factor,DamagePenaltyFactorEnum);
     803        parameters->FindParam(&max_damage,DamageMaxDamageEnum);
     804
     805        //right hand side penalizes to max_damage
     806        pe->values[0]=kmax*pow(10.,penalty_factor)*max_damage;
     807
     808        /*Clean up and return*/
     809        return pe;
     810}
     811/*}}}*/
     812#endif
     813
    697814/*FUNCTION Pengrid::ResetConstraint {{{*/
    698815void  Pengrid::ResetConstraint(void){

  • issm/trunk-jpl/src/c/classes/Loads/Pengrid.h

    r16125 r16188  
    9191                void           ConstraintActivateThermal(int* punstable);
    9292                #endif
     93                #ifdef _HAVE_DAMAGE_
     94                ElementMatrix* PenaltyCreateKMatrixDamageEvolution(IssmDouble kmax);
     95                ElementVector* PenaltyCreatePVectorDamageEvolution(IssmDouble kmax);
     96                void           ConstraintActivateDamageEvolution(int* punstable);
     97                #endif
     98
    9399                #ifdef _HAVE_HYDROLOGY_
    94100                ElementMatrix* PenaltyCreateKMatrixHydrologyDCInefficient(IssmDouble kmax);

  • issm/trunk-jpl/src/c/modules/ConstraintsStatex/ConstraintsStateLocal.h

    r15000 r16188  
    66#define _CONSTRAINTSSTATELOCAL_H
    77
     8
     9#ifdef HAVE_CONFIG_H
     10        #include <config.h>
     11#else
     12#error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!"
     13#endif
     14
    815#include "../../classes/classes.h"
    916
    1017/*melting: */
    1118void  ThermalConstraintsState(Loads* loads, int* pconverged, int* pnum_unstable_constraints,int analysis_type);
     19
     20/*damage: */
     21#ifdef _HAVE_DAMAGE_
     22void  DamageConstraintsState(Loads* loads, int* pconverged, int* pnum_unstable_constraints,int analysis_type);
     23#endif
    1224
    1325/*rifts module: */

  • issm/trunk-jpl/src/c/modules/ConstraintsStatex/ConstraintsStatex.cpp

    r15849 r16188  
    3232        #ifdef _HAVE_RIFTS_
    3333        if(RiftIsPresent(femmodel->loads,analysis_type)){
     34
    3435                RiftConstraintsState(&converged,&num_unstable_constraints,femmodel->loads,min_mechanical_constraints,analysis_type);
    3536        }
    3637        #endif
     38
    3739        #ifdef _HAVE_THERMAL_
    38         if(ThermalIsPresent(femmodel->loads,analysis_type)){
     40
     41        if(PengridIsPresent(femmodel->loads,ThermalAnalysisEnum)){
     42
    3943                ThermalConstraintsState(femmodel->loads,&converged,&num_unstable_constraints,analysis_type);
    4044        }
    4145        #endif
     46
     47        #ifdef _HAVE_DAMAGE_
     48
     49        if(PengridIsPresent(femmodel->loads,DamageEvolutionAnalysisEnum)){
     50
     51                DamageConstraintsState(femmodel->loads,&converged,&num_unstable_constraints,analysis_type);
     52        }
     53        #endif
     54
    4255
    4356        /*Assign output pointers: */

  • issm/trunk-jpl/src/c/modules/ConstraintsStatex/ConstraintsStatex.h

    r15849 r16188  
    99
    1010/* local prototypes: */
    11 int  ThermalIsPresent(Loads* loads,int analysis_type);
     11int  PengridIsPresent(Loads* loads,int analysis_type);
    1212int  RiftIsPresent(Loads* loads,int analysis_type);
    1313void ConstraintsStatex(int* pconverged, int* pnum_unstable_constraints,FemModel* femmodel);

  • issm/trunk-jpl/src/c/modules/ModelProcessorx/Damage/CreateLoadsDamage.cpp

    r16181 r16188  
    88
    99void    CreateLoadsDamage(Loads** ploads, IoModel* iomodel){
     10
     11        /*Recover pointer: */
     12        Loads* loads=*ploads;
     13
     14        /*create penalties for nodes: no node can have a damage > 1*/
     15        iomodel->FetchData(1,DamageSpcdamageEnum);
     16        CreateSingleNodeToElementConnectivity(iomodel);
     17
     18        for(int i=0;i<iomodel->numberofvertices;i++){
     19
     20                /*keep only this partition's nodes:*/
     21                if((iomodel->my_vertices[i]==1)){
     22                        if (xIsNan<IssmDouble>(iomodel->Data(DamageSpcdamageEnum)[i])){ //No penalty applied on spc nodes!
     23                                loads->AddObject(new Pengrid(iomodel->loadcounter+i+1,i,iomodel,DamageEvolutionAnalysisEnum));
     24                        }
     25                }
     26        }
     27        iomodel->DeleteData(1,DamageSpcdamageEnum);
     28
     29        /*Assign output pointer: */
     30        *ploads=loads;
    1031       
    11         /*No loads*/
    12 
    1332}

  • issm/trunk-jpl/src/c/modules/ModelProcessorx/Damage/CreateParametersDamage.cpp

    r16181 r16188  
    2121       
    2222        iomodel->Constant(&law,DamageLawEnum);
     23               
     24        parameters->AddObject(iomodel->CopyConstantObject(DamagePenaltyThresholdEnum));
     25        parameters->AddObject(iomodel->CopyConstantObject(DamagePenaltyLockEnum));
     26        parameters->AddObject(iomodel->CopyConstantObject(DamagePenaltyFactorEnum));
     27        parameters->AddObject(iomodel->CopyConstantObject(DamageMaxiterEnum));
     28        parameters->AddObject(iomodel->CopyConstantObject(DamageMaxDamageEnum));
    2329
    2430        /*Retrieve law dependent parameters: */

  • issm/trunk-jpl/src/c/modules/ModelProcessorx/Damage/UpdateElementsDamage.cpp

    r16181 r16188  
    2626        iomodel->FetchDataToInput(elements,VzEnum);
    2727        iomodel->FetchDataToInput(elements,DamageDEnum);
     28        iomodel->FetchDataToInput(elements,MaskIceLevelsetEnum);
     29        iomodel->FetchDataToInput(elements,PressureEnum);
    2830
    2931        bool dakota_analysis;

  • issm/trunk-jpl/src/c/modules/ResetConstraintsx/ResetConstraintsx.cpp

    r16144 r16188  
    3131        #endif
    3232        #ifdef _HAVE_THERMAL_
    33         if(ThermalIsPresent(femmodel->loads,analysis_type)){
     33        if(PengridIsPresent(femmodel->loads,ThermalAnalysisEnum)){
    3434                ThermalConstraintsReset(femmodel->loads,analysis_type);
    3535        }
    3636        #endif
     37        #ifdef _HAVE_DAMAGE_
     38        if(PengridIsPresent(femmodel->loads,DamageEvolutionAnalysisEnum)){
     39                DamageConstraintsReset(femmodel->loads,analysis_type);
     40        }
     41        #endif
     42
    3743}

  • issm/trunk-jpl/src/c/modules/ResetConstraintsx/ResetConstraintsx.h

    r15849 r16188  
    1010/* local prototypes: */
    1111void ThermalConstraintsReset(Loads* loads, int configuration_type);
     12void DamageConstraintsReset(Loads* loads, int configuration_type);
    1213void ResetConstraintsx(FemModel* femmodel);
    1314

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

    r16181 r16188  
    163163        DamageStressThresholdEnum,
    164164        DamageStabilizationEnum,
     165        DamagePenaltyThresholdEnum,
     166        DamagePenaltyLockEnum,
     167        DamagePenaltyFactorEnum,
     168        DamageMaxiterEnum,
    165169        DamageSpcdamageEnum,
     170        DamageMaxDamageEnum,
    166171        MaterialsRhoIceEnum,
    167172        MaterialsRhoWaterEnum,

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

    r16181 r16188  
    171171                case DamageStressThresholdEnum : return "DamageStressThreshold";
    172172                case DamageStabilizationEnum : return "DamageStabilization";
     173                case DamagePenaltyThresholdEnum : return "DamagePenaltyThreshold";
     174                case DamagePenaltyLockEnum : return "DamagePenaltyLock";
     175                case DamagePenaltyFactorEnum : return "DamagePenaltyFactor";
     176                case DamageMaxiterEnum : return "DamageMaxiter";
    173177                case DamageSpcdamageEnum : return "DamageSpcdamage";
     178                case DamageMaxDamageEnum : return "DamageMaxDamage";
    174179                case MaterialsRhoIceEnum : return "MaterialsRhoIce";
    175180                case MaterialsRhoWaterEnum : return "MaterialsRhoWater";

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

    r16181 r16188  
    174174              else if (strcmp(name,"DamageStressThreshold")==0) return DamageStressThresholdEnum;
    175175              else if (strcmp(name,"DamageStabilization")==0) return DamageStabilizationEnum;
     176              else if (strcmp(name,"DamagePenaltyThreshold")==0) return DamagePenaltyThresholdEnum;
     177              else if (strcmp(name,"DamagePenaltyLock")==0) return DamagePenaltyLockEnum;
     178              else if (strcmp(name,"DamagePenaltyFactor")==0) return DamagePenaltyFactorEnum;
     179              else if (strcmp(name,"DamageMaxiter")==0) return DamageMaxiterEnum;
    176180              else if (strcmp(name,"DamageSpcdamage")==0) return DamageSpcdamageEnum;
     181              else if (strcmp(name,"DamageMaxDamage")==0) return DamageMaxDamageEnum;
    177182              else if (strcmp(name,"MaterialsRhoIce")==0) return MaterialsRhoIceEnum;
    178183              else if (strcmp(name,"MaterialsRhoWater")==0) return MaterialsRhoWaterEnum;
    … …  
    255260              else if (strcmp(name,"SurfaceforcingsHref")==0) return SurfaceforcingsHrefEnum;
    256261              else if (strcmp(name,"SurfaceforcingsSmbref")==0) return SurfaceforcingsSmbrefEnum;
    257               else if (strcmp(name,"SurfaceforcingsBPos")==0) return SurfaceforcingsBPosEnum;
     262         else stage=3;
     263   }
     264   if(stage==3){
     265              if (strcmp(name,"SurfaceforcingsBPos")==0) return SurfaceforcingsBPosEnum;
    258266              else if (strcmp(name,"SurfaceforcingsBNeg")==0) return SurfaceforcingsBNegEnum;
    259267              else if (strcmp(name,"ThermalIsenthalpy")==0) return ThermalIsenthalpyEnum;
    260268              else if (strcmp(name,"ThermalIsdynamicbasalspc")==0) return ThermalIsdynamicbasalspcEnum;
    261269              else if (strcmp(name,"ThermalMaxiter")==0) return ThermalMaxiterEnum;
    262          else stage=3;
    263    }
    264    if(stage==3){
    265               if (strcmp(name,"ThermalPenaltyFactor")==0) return ThermalPenaltyFactorEnum;
     270              else if (strcmp(name,"ThermalPenaltyFactor")==0) return ThermalPenaltyFactorEnum;
    266271              else if (strcmp(name,"ThermalPenaltyLock")==0) return ThermalPenaltyLockEnum;
    267272              else if (strcmp(name,"ThermalPenaltyThreshold")==0) return ThermalPenaltyThresholdEnum;
    … …  
    378383              else if (strcmp(name,"IntParam")==0) return IntParamEnum;
    379384              else if (strcmp(name,"IntVecParam")==0) return IntVecParamEnum;
    380               else if (strcmp(name,"TransientParam")==0) return TransientParamEnum;
     385         else stage=4;
     386   }
     387   if(stage==4){
     388              if (strcmp(name,"TransientParam")==0) return TransientParamEnum;
    381389              else if (strcmp(name,"IceFrontType")==0) return IceFrontTypeEnum;
    382390              else if (strcmp(name,"SSA2dIceFront")==0) return SSA2dIceFrontEnum;
    383391              else if (strcmp(name,"SSA3dIceFront")==0) return SSA3dIceFrontEnum;
    384392              else if (strcmp(name,"Matice")==0) return MaticeEnum;
    385          else stage=4;
    386    }
    387    if(stage==4){
    388               if (strcmp(name,"Matpar")==0) return MatparEnum;
     393              else if (strcmp(name,"Matpar")==0) return MatparEnum;
    389394              else if (strcmp(name,"Node")==0) return NodeEnum;
    390395              else if (strcmp(name,"Numericalflux")==0) return NumericalfluxEnum;
    … …  
    501506              else if (strcmp(name,"BasalFriction")==0) return BasalFrictionEnum;
    502507              else if (strcmp(name,"ViscousHeating")==0) return ViscousHeatingEnum;
    503               else if (strcmp(name,"QmuTemperature")==0) return QmuTemperatureEnum;
     508         else stage=5;
     509   }
     510   if(stage==5){
     511              if (strcmp(name,"QmuTemperature")==0) return QmuTemperatureEnum;
    504512              else if (strcmp(name,"HydrologyWaterVx")==0) return HydrologyWaterVxEnum;
    505513              else if (strcmp(name,"HydrologyWaterVy")==0) return HydrologyWaterVyEnum;
    506514              else if (strcmp(name,"StressTensor")==0) return StressTensorEnum;
    507515              else if (strcmp(name,"StressTensorxx")==0) return StressTensorxxEnum;
    508          else stage=5;
    509    }
    510    if(stage==5){
    511               if (strcmp(name,"StressTensorxy")==0) return StressTensorxyEnum;
     516              else if (strcmp(name,"StressTensorxy")==0) return StressTensorxyEnum;
    512517              else if (strcmp(name,"StressTensorxz")==0) return StressTensorxzEnum;
    513518              else if (strcmp(name,"StressTensoryy")==0) return StressTensoryyEnum;

  • issm/trunk-jpl/src/c/solutionsequences/solutionsequences.h

    r15564 r16188  
    1313
    1414void solutionsequence_thermal_nonlinear(FemModel* femmodel);
     15void solutionsequence_damage_nonlinear(FemModel* femmodel);
    1516void solutionsequence_hydro_nonlinear(FemModel* femmodel);
    1617void solutionsequence_nonlinear(FemModel* femmodel,bool conserve_loads);

  • issm/trunk-jpl/src/m/classes/damage.m

    r16181 r16188  
    1111                law                 = '';
    1212                spcdamage           = NaN;
     13                max_damage          = NaN;
     14       
     15                %numerical
    1316                stabilization       = NaN;
     17                penalty_threshold   = NaN;
     18                maxiter             = NaN;
     19                penalty_lock        = NaN;
     20                penalty_factor      = NaN;
    1421               
    1522                %parameters for law 'initial':
    … …  
    4653                        obj.D=0;
    4754                        obj.law='undamaged';
     55                       
     56                        obj.max_damage=1-1e-5; %if damage reaches 1, solve becomes singular, as viscosity becomes nil
    4857               
    4958                        %Type of stabilization used
    5059                        obj.stabilization=2;
    51        
     60                       
     61                        %Maximum number of iterations
     62                        obj.maxiter=100;
     63
     64                        %factor used to compute the values of the penalties: kappa=max(stiffness matrix)*10^penalty_factor
     65                        obj.penalty_factor=3;
     66                       
     67                        %stabilize unstable damage constraints that keep zigzagging after n iteration (default is 0, no stabilization)
     68                        obj.penalty_lock=0;
     69                       
     70                        %threshold to declare convergence of thermal solution (default is 0)
     71                        obj.penalty_threshold=0;
     72               
    5273                        %damage parameters for 'initial' law of damage propagation
    5374                        obj.stress_threshold=0;
    … …  
    6182                function md = checkconsistency(obj,md,solution,analyses) % {{{
    6283                       
    63                         md = checkfield(md,'damage.D','>=0',0,'size',[md.mesh.numberofvertices 1]);
     84                        md = checkfield(md,'damage.D','>=',0,'<=',obj.max_damage,'size',[md.mesh.numberofvertices 1]);
     85                        md = checkfield(md,'damage.max_damage','<',1,'>=',0);
    6486                        md = checkfield(md,'damage.law','values',{'undamaged','pralong'});
    6587                        md = checkfield(md,'damage.spcdamage','forcing',1);
     88                       
    6689                        md = checkfield(md,'damage.stabilization','numel',[1],'values',[0 1 2]);
     90                        md = checkfield(md,'damage.maxiter','>=0',0);
     91                        md = checkfield(md,'damage.penalty_factor','>=0',0);
     92                        md = checkfield(md,'damage.penalty_lock','>=0',0);
     93                        md = checkfield(md,'damage.penalty_threshold','>=0',0);
     94
    6795                        if strcmpi(obj.law,'pralong'),
    6896                                md = checkfield(md,'damage.healing','>=',0);
    … …  
    72100                                md = checkfield(md,'damage.c4','>=',0);
    73101                                md = checkfield(md,'damage.stress_threshold','>=',0);
     102                        elseif strcmpi(obj.law,'undamaged'),
     103                                if (solution==DamageEvolutionSolutionEnum),
     104                                        error('you are trying to evolve an undamaged material. Choose a valid evolution law (md.damage.law)');
     105                                end
     106                        else
     107                                error('invalid damage evolution law');
    74108                        end
    75109
    … …  
    81115                        fielddisplay(obj,'law','damage law (string) from {''undamaged'',''pralong''}');
    82116                        fielddisplay(obj,'spcdamage','damage constraints (NaN means no constraint) [m]');
     117                        fielddisplay(obj,'max_damage','maximum damage sustained by ice (0<=max_damage<1)');
     118                       
    83119                        fielddisplay(obj,'stabilization','0: no, 1: artificial_diffusivity, 2: SUPG');
     120                        fielddisplay(obj,'maxiter','maximum number of non linear iterations');
     121                        fielddisplay(obj,'penalty_lock','stabilize unstable damage constraints that keep zigzagging after n iteration (default is 0, no stabilization)');
     122                        fielddisplay(obj,'penalty_threshold','threshold to declare convergence of damage evolution solution (default is 0)');
     123                        fielddisplay(obj,'penalty_factor','scaling exponent (default is 3)');
     124
    84125                        if strcmpi(obj.law,'pralong'),
    85126                                fielddisplay(obj,'c1','damage parameter 1');
    … …  
    96137                        WriteData(fid,'object',obj,'fieldname','D','format','DoubleMat','mattype',1);
    97138                        WriteData(fid,'object',obj,'fieldname','law','format','String');
     139                        WriteData(fid,'object',obj,'fieldname','spcdamage','format','DoubleMat','mattype',1,'forcinglength',md.mesh.numberofvertices+1);
     140                        WriteData(fid,'object',obj,'fieldname','max_damage','format','Double');
     141
    98142                        WriteData(fid,'object',obj,'fieldname','stabilization','format','Integer');
    99                         WriteData(fid,'object',obj,'fieldname','spcdamage','format','DoubleMat','mattype',1,'forcinglength',md.mesh.numberofvertices+1);
     143                        WriteData(fid,'object',obj,'fieldname','penalty_threshold','format','Integer');
     144                        WriteData(fid,'object',obj,'fieldname','maxiter','format','Integer');
     145                        WriteData(fid,'object',obj,'fieldname','penalty_lock','format','Integer');
     146                        WriteData(fid,'object',obj,'fieldname','penalty_factor','format','Double');
     147       
    100148                        if strcmpi(obj.law,'pralong'),
    101149                                WriteData(fid,'object',obj,'fieldname','c1','format','Double');

  • issm/trunk-jpl/src/m/classes/thermal.m

    r16144 r16188  
    6767                        fielddisplay(obj,'penalty_lock','stabilize unstable thermal constraints that keep zigzagging after n iteration (default is 0, no stabilization)');
    6868                        fielddisplay(obj,'penalty_threshold','threshold to declare convergence of thermal solution (default is 0)');
     69                        fielddisplay(obj,'penalty_factor','scaling exponent (default is 3)');
    6970                        fielddisplay(obj,'isenthalpy','use an enthalpy formulation to include temperate ice (default is 0)');
    7071                        fielddisplay(obj,'isdynamicbasalspc',['enable dynamic setting of basal forcing. required for enthalpy formulation (default is 0)']);

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

    r16181 r16188  
    163163def DamageStressThresholdEnum(): return StringToEnum("DamageStressThreshold")[0]
    164164def DamageStabilizationEnum(): return StringToEnum("DamageStabilization")[0]
     165def DamagePenaltyThresholdEnum(): return StringToEnum("DamagePenaltyThreshold")[0]
     166def DamagePenaltyLockEnum(): return StringToEnum("DamagePenaltyLock")[0]
     167def DamagePenaltyFactorEnum(): return StringToEnum("DamagePenaltyFactor")[0]
     168def DamageMaxiterEnum(): return StringToEnum("DamageMaxiter")[0]
    165169def DamageSpcdamageEnum(): return StringToEnum("DamageSpcdamage")[0]
     170def DamageMaxDamageEnum(): return StringToEnum("DamageMaxDamage")[0]
    166171def MaterialsRhoIceEnum(): return StringToEnum("MaterialsRhoIce")[0]
    167172def MaterialsRhoWaterEnum(): return StringToEnum("MaterialsRhoWater")[0]

  • issm/trunk-jpl/test/NightlyRun/test271.m

    r16181 r16188  
    44md.damage.D=0.5*ones(md.mesh.numberofvertices,1);
    55md.damage.spcdamage=NaN*ones(md.mesh.numberofvertices,1);
     6md.damage.law='pralong';
    67
    78%boundary conditions for damage, to be put in SquareShelf.par
    … …  
    1112md.damage.spcdamage(pos)=0;
    1213
    13 md.verbose=verbose('11111111');
     14md.verbose=verbose('solution',true);
    1415
    1516md=setflowequation(md,'SSA','all');
Note: See TracChangeset for help on using the changeset viewer.