Changeset 18194


Ignore:
Timestamp:
06/27/14 14:42:12 (11 years ago)
Author:
cborstad
Message:

NEW: working on damage calculation from stress balance solution

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

Legend:

Unmodified
Added
Removed

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

    r18193 r18194  
    4141        this->inputs->AddInput(input_in);
    4242}/*}}}*/
     43void       Element::ComputeNewDamage(){/*{{{*/
     44
     45        IssmDouble *xyz_list=NULL;
     46        IssmDouble  eps_xx,eps_xy,eps_yy,eps_xz,eps_yz,eps_zz,eps_eff;
     47        IssmDouble  epsmin=1.e-27;
     48        int         dim;
     49
     50        /*Retrieve all inputs we will be needing: */
     51        /* TODO: retrieve parameters eps_0 and eps_f and input DamageD(bar?) */
     52        this->GetVerticesCoordinates(&xyz_list);
     53        this->ComputeStrainRate();
     54        parameters->FindParam(&dim,DomainDimensionEnum);
     55        Input* eps_xx_input=this->GetInput(StrainRatexxEnum); _assert_(eps_xx_input);
     56        Input* eps_yy_input=this->GetInput(StrainRateyyEnum); _assert_(eps_yy_input);
     57        Input* eps_xy_input=this->GetInput(StrainRatexyEnum); _assert_(eps_xy_input);
     58        Input* eps_xz_input=NULL;
     59        Input* eps_yz_input=NULL;
     60        Input* eps_zz_input=NULL;
     61        if(dim==3){
     62                eps_xz_input=this->GetInput(StrainRatexzEnum); _assert_(eps_xz_input);
     63                eps_yz_input=this->GetInput(StrainRateyzEnum); _assert_(eps_yz_input);
     64                eps_zz_input=this->GetInput(StrainRatezzEnum); _assert_(eps_zz_input);
     65        }
     66
     67        /* Start looping on the number of vertices: */
     68        Gauss* gauss=this->NewGauss();
     69        int numvertices = this->GetNumberOfVertices();
     70        for (int iv=0;iv<numvertices;iv++){
     71                gauss->GaussVertex(iv);
     72
     73                eps_xx_input->GetInputValue(&eps_xx,gauss);
     74                eps_yy_input->GetInputValue(&eps_yy,gauss);
     75                eps_xy_input->GetInputValue(&eps_xy,gauss);
     76                if(dim==3){
     77                        eps_xz_input->GetInputValue(&eps_xz,gauss);
     78                        eps_yz_input->GetInputValue(&eps_yz,gauss);
     79                        eps_zz_input->GetInputValue(&eps_zz,gauss);
     80                }
     81                else{eps_xz=0; eps_yz=0; eps_zz=0;}
     82
     83                /* eps_eff^2 = exx^2 + eyy^2 + exy^2 + exz^2 + eyz^2 + exx*eyy */
     84                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);
     85
     86                /*TODO: compute kappa from initial D, then compute new D */
     87
     88        }
     89
     90        /* TODO: add newdamage input to DamageEnum and NewDamageEnum */
     91
     92        /*Clean up and return*/
     93        xDelete<IssmDouble>(xyz_list);
     94        delete gauss;
     95
     96}/*}}}*/
     97void       Element::ComputeStrainRate(){/*{{{*/
     98
     99        int         dim;
     100        IssmDouble *xyz_list = NULL;
     101        IssmDouble  epsilon[6];
     102
     103        /*Retrieve all inputs we will be needing: */
     104        this->GetVerticesCoordinates(&xyz_list);
     105        parameters->FindParam(&dim,DomainDimensionEnum);
     106        Input* vx_input=this->GetInput(VxEnum); _assert_(vx_input);
     107        Input* vy_input=this->GetInput(VyEnum); _assert_(vy_input);
     108        Input* vz_input=NULL;
     109        if(dim==3){vz_input=this->GetInput(VzEnum); _assert_(vz_input);}
     110
     111        /*Allocate arrays*/
     112        int numvertices = this->GetNumberOfVertices();
     113        IssmDouble* eps_xx = xNew<IssmDouble>(numvertices);
     114        IssmDouble* eps_yy = xNew<IssmDouble>(numvertices);
     115        IssmDouble* eps_zz = xNew<IssmDouble>(numvertices);
     116        IssmDouble* eps_xy = xNew<IssmDouble>(numvertices);
     117        IssmDouble* eps_xz = xNew<IssmDouble>(numvertices);
     118        IssmDouble* eps_yz = xNew<IssmDouble>(numvertices);
     119
     120        /* Start looping on the number of vertices: */
     121        Gauss* gauss=this->NewGauss();
     122        for (int iv=0;iv<numvertices;iv++){
     123                gauss->GaussVertex(iv);
     124
     125                /*Compute strain rate viscosity and pressure: */
     126                if(dim==2)
     127                 this->StrainRateSSA(&epsilon[0],xyz_list,gauss,vx_input,vy_input);
     128                else
     129                 this->StrainRateFS(&epsilon[0],xyz_list,gauss,vx_input,vy_input,vz_input);
     130
     131                if(dim==2){
     132                         /* epsilon=[exx,eyy,exy];*/
     133                        eps_xx[iv]=epsilon[0];
     134                        eps_yy[iv]=epsilon[1];
     135                        eps_xy[iv]=epsilon[2];
     136                }
     137                else{
     138                        /*epsilon=[exx eyy ezz exy exz eyz]*/
     139                        eps_xx[iv]=epsilon[0];
     140                        eps_yy[iv]=epsilon[1];
     141                        eps_zz[iv]=epsilon[2];
     142                        eps_xy[iv]=epsilon[3];
     143                        eps_xz[iv]=epsilon[4];
     144                        eps_yz[iv]=epsilon[5];
     145                }
     146        }
     147
     148        /*Add Stress tensor components into inputs*/
     149        this->AddInput(StrainRatexxEnum,eps_xx,P1Enum);
     150        this->AddInput(StrainRatexyEnum,eps_xy,P1Enum);
     151        this->AddInput(StrainRatexzEnum,eps_xz,P1Enum);
     152        this->AddInput(StrainRateyyEnum,eps_yy,P1Enum);
     153        this->AddInput(StrainRateyzEnum,eps_yz,P1Enum);
     154        this->AddInput(StrainRatezzEnum,eps_zz,P1Enum);
     155
     156        /*Clean up and return*/
     157        delete gauss;
     158        xDelete<IssmDouble>(xyz_list);
     159        xDelete<IssmDouble>(eps_xx);
     160        xDelete<IssmDouble>(eps_yy);
     161        xDelete<IssmDouble>(eps_zz);
     162        xDelete<IssmDouble>(eps_xy);
     163        xDelete<IssmDouble>(eps_xz);
     164        xDelete<IssmDouble>(eps_yz);
     165
     166}
     167/*}}}*/
    43168void       Element::CoordinateSystemTransform(IssmDouble** ptransform,Node** nodes_list,int numnodes,int* cs_array){/*{{{*/
    44169
    … …  
    113238        /*Assign output pointer*/
    114239        *ptransform=transform;
    115 }
    116 /*}}}*/
    117 void       Element::ComputeStrainRate(){/*{{{*/
    118 
    119         int         dim;
    120         IssmDouble *xyz_list = NULL;
    121         IssmDouble  epsilon[6];
    122 
    123         /*Retrieve all inputs we will be needing: */
    124         this->GetVerticesCoordinates(&xyz_list);
    125         parameters->FindParam(&dim,DomainDimensionEnum);
    126         Input* vx_input=this->GetInput(VxEnum); _assert_(vx_input);
    127         Input* vy_input=this->GetInput(VyEnum); _assert_(vy_input);
    128         Input* vz_input=NULL;
    129         if(dim==3){vz_input=this->GetInput(VzEnum); _assert_(vz_input);}
    130 
    131         /*Allocate arrays*/
    132         int numvertices = this->GetNumberOfVertices();
    133         IssmDouble* eps_xx = xNew<IssmDouble>(numvertices);
    134         IssmDouble* eps_yy = xNew<IssmDouble>(numvertices);
    135         IssmDouble* eps_zz = xNew<IssmDouble>(numvertices);
    136         IssmDouble* eps_xy = xNew<IssmDouble>(numvertices);
    137         IssmDouble* eps_xz = xNew<IssmDouble>(numvertices);
    138         IssmDouble* eps_yz = xNew<IssmDouble>(numvertices);
    139 
    140         /* Start looping on the number of vertices: */
    141         Gauss* gauss=this->NewGauss();
    142         for (int iv=0;iv<numvertices;iv++){
    143                 gauss->GaussVertex(iv);
    144 
    145                 /*Compute strain rate viscosity and pressure: */
    146                 if(dim==2)
    147                  this->StrainRateSSA(&epsilon[0],xyz_list,gauss,vx_input,vy_input);
    148                 else
    149                  this->StrainRateFS(&epsilon[0],xyz_list,gauss,vx_input,vy_input,vz_input);
    150 
    151                 if(dim==2){
    152                          /* epsilon=[exx,eyy,exy];*/
    153                         eps_xx[iv]=epsilon[0];
    154                         eps_yy[iv]=epsilon[1];
    155                         eps_xy[iv]=epsilon[2];
    156                 }
    157                 else{
    158                         /*epsilon=[exx eyy ezz exy exz eyz]*/
    159                         eps_xx[iv]=epsilon[0];
    160                         eps_yy[iv]=epsilon[1];
    161                         eps_zz[iv]=epsilon[2];
    162                         eps_xy[iv]=epsilon[3];
    163                         eps_xz[iv]=epsilon[4];
    164                         eps_yz[iv]=epsilon[5];
    165                 }
    166         }
    167 
    168         /*Add Stress tensor components into inputs*/
    169         this->AddInput(StrainRatexxEnum,eps_xx,P1Enum);
    170         this->AddInput(StrainRatexyEnum,eps_xy,P1Enum);
    171         this->AddInput(StrainRatexzEnum,eps_xz,P1Enum);
    172         this->AddInput(StrainRateyyEnum,eps_yy,P1Enum);
    173         this->AddInput(StrainRateyzEnum,eps_yz,P1Enum);
    174         this->AddInput(StrainRatezzEnum,eps_zz,P1Enum);
    175 
    176         /*Clean up and return*/
    177         delete gauss;
    178         xDelete<IssmDouble>(xyz_list);
    179         xDelete<IssmDouble>(eps_xx);
    180         xDelete<IssmDouble>(eps_yy);
    181         xDelete<IssmDouble>(eps_zz);
    182         xDelete<IssmDouble>(eps_xy);
    183         xDelete<IssmDouble>(eps_xz);
    184         xDelete<IssmDouble>(eps_yz);
    185 
    186240}
    187241/*}}}*/
    … …  
    10531107                                input=this->inputs->GetInput(output_enum);
    10541108                                break;
     1109                        case NewDamageEnum:
     1110                                this->ComputeNewDamage();
     1111                                input=this->inputs->GetInput(output_enum);
     1112                                break;
    10551113                        default:
    10561114                                _error_("input "<<EnumToStringx(output_enum)<<" not found in element");

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

    r18192 r18194  
    5959                /* bool       AllActive(void); */
    6060                /* bool       AnyActive(void); */
     61                void       ComputeNewDamage();
     62                void       ComputeStrainRate();
    6163                void       CoordinateSystemTransform(IssmDouble** ptransform,Node** nodes,int numnodes,int* cs_array);
    62                 void       ComputeStrainRate();
    6364                void       Echo();
    6465                void       DeepEcho();

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

    r18179 r18194  
    200200        DamageEvolutionNumRequestedOutputsEnum,
    201201        DamageEvolutionRequestedOutputsEnum,
     202        NewDamageEnum,
    202203        MaterialsRhoIceEnum,
    203204        MaterialsRhoSeawaterEnum,

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

    r18179 r18194  
    208208                case DamageEvolutionNumRequestedOutputsEnum : return "DamageEvolutionNumRequestedOutputs";
    209209                case DamageEvolutionRequestedOutputsEnum : return "DamageEvolutionRequestedOutputs";
     210                case NewDamageEnum : return "NewDamage";
    210211                case MaterialsRhoIceEnum : return "MaterialsRhoIce";
    211212                case MaterialsRhoSeawaterEnum : return "MaterialsRhoSeawater";

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

    r18179 r18194  
    211211              else if (strcmp(name,"DamageEvolutionNumRequestedOutputs")==0) return DamageEvolutionNumRequestedOutputsEnum;
    212212              else if (strcmp(name,"DamageEvolutionRequestedOutputs")==0) return DamageEvolutionRequestedOutputsEnum;
     213              else if (strcmp(name,"NewDamage")==0) return NewDamageEnum;
    213214              else if (strcmp(name,"MaterialsRhoIce")==0) return MaterialsRhoIceEnum;
    214215              else if (strcmp(name,"MaterialsRhoSeawater")==0) return MaterialsRhoSeawaterEnum;
    … …  
    259260              else if (strcmp(name,"MassFluxSegmentsPresent")==0) return MassFluxSegmentsPresentEnum;
    260261              else if (strcmp(name,"QmuMassFluxSegmentsPresent")==0) return QmuMassFluxSegmentsPresentEnum;
    261               else if (strcmp(name,"QmuNumberofpartitions")==0) return QmuNumberofpartitionsEnum;
    262262         else stage=3;
    263263   }
    264264   if(stage==3){
    265               if (strcmp(name,"QmuNumberofresponses")==0) return QmuNumberofresponsesEnum;
     265              if (strcmp(name,"QmuNumberofpartitions")==0) return QmuNumberofpartitionsEnum;
     266              else if (strcmp(name,"QmuNumberofresponses")==0) return QmuNumberofresponsesEnum;
    266267              else if (strcmp(name,"QmuPartition")==0) return QmuPartitionEnum;
    267268              else if (strcmp(name,"QmuResponsedescriptors")==0) return QmuResponsedescriptorsEnum;
    … …  
    382383              else if (strcmp(name,"HydrologyDCInefficientAnalysis")==0) return HydrologyDCInefficientAnalysisEnum;
    383384              else if (strcmp(name,"HydrologyDCEfficientAnalysis")==0) return HydrologyDCEfficientAnalysisEnum;
    384               else if (strcmp(name,"HydrologySolution")==0) return HydrologySolutionEnum;
    385385         else stage=4;
    386386   }
    387387   if(stage==4){
    388               if (strcmp(name,"MeltingAnalysis")==0) return MeltingAnalysisEnum;
     388              if (strcmp(name,"HydrologySolution")==0) return HydrologySolutionEnum;
     389              else if (strcmp(name,"MeltingAnalysis")==0) return MeltingAnalysisEnum;
    389390              else if (strcmp(name,"MasstransportAnalysis")==0) return MasstransportAnalysisEnum;
    390391              else if (strcmp(name,"MasstransportSolution")==0) return MasstransportSolutionEnum;
    … …  
    505506              else if (strcmp(name,"Fill")==0) return FillEnum;
    506507              else if (strcmp(name,"FractionIncrement")==0) return FractionIncrementEnum;
    507               else if (strcmp(name,"Friction")==0) return FrictionEnum;
    508508         else stage=5;
    509509   }
    510510   if(stage==5){
    511               if (strcmp(name,"Internal")==0) return InternalEnum;
     511              if (strcmp(name,"Friction")==0) return FrictionEnum;
     512              else if (strcmp(name,"Internal")==0) return InternalEnum;
    512513              else if (strcmp(name,"MassFlux")==0) return MassFluxEnum;
    513514              else if (strcmp(name,"MeltingOffset")==0) return MeltingOffsetEnum;
    … …  
    628629              else if (strcmp(name,"Step")==0) return StepEnum;
    629630              else if (strcmp(name,"Time")==0) return TimeEnum;
    630               else if (strcmp(name,"WaterColumnOld")==0) return WaterColumnOldEnum;
    631631         else stage=6;
    632632   }
    633633   if(stage==6){
    634               if (strcmp(name,"Outputdefinition")==0) return OutputdefinitionEnum;
     634              if (strcmp(name,"WaterColumnOld")==0) return WaterColumnOldEnum;
     635              else if (strcmp(name,"Outputdefinition")==0) return OutputdefinitionEnum;
    635636              else if (strcmp(name,"OutputdefinitionList")==0) return OutputdefinitionListEnum;
    636637              else if (strcmp(name,"Massfluxatgate")==0) return MassfluxatgateEnum;

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

    r18179 r18194  
    200200def DamageEvolutionNumRequestedOutputsEnum(): return StringToEnum("DamageEvolutionNumRequestedOutputs")[0]
    201201def DamageEvolutionRequestedOutputsEnum(): return StringToEnum("DamageEvolutionRequestedOutputs")[0]
     202def NewDamageEnum(): return StringToEnum("NewDamage")[0]
    202203def MaterialsRhoIceEnum(): return StringToEnum("MaterialsRhoIce")[0]
    203204def MaterialsRhoSeawaterEnum(): return StringToEnum("MaterialsRhoSeawater")[0]

  • TabularUnified issm/trunk-jpl/src/m/plot/plotmodel.py

    r17919 r18194  
    5454        if numberofplots:
    5555               
    56                 if plt.fignum_exists(figurenumber):
    57                         plt.cla()
     56                #if plt.fignum_exists(figurenumber):
     57                #       plt.cla()
    5858
    5959                #if figsize specified
Note: See TracChangeset for help on using the changeset viewer.