Changeset 17757


Ignore:
Timestamp:
04/17/14 13:15:15 (11 years ago)
Author:
cborstad
Message:

CHG: stripping damage out of solution if MaterialsType is MaticeEnum

Location:
issm/trunk-jpl/src
Files:
2 added
13 edited

Legend:

Unmodified
Added
Removed

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

    r17700 r17757  
    131131        bool   dakota_analysis;
    132132        bool   islevelset;
     133        bool   isdamage;
    133134
    134135        /*Fetch constants needed: */
    … …  
    154155        else
    155156         iscoupling = false;
     157
     158        /*is damage mechanics being used?*/
     159        if(materials_type==MaticeEnum) isdamage = false;
     160        else if(materials_type==MatdamageiceEnum) isdamage = true;
     161        else _error_("Material type not recognized");
    156162
    157163        /*Get finite element type*/
    … …  
    212218        iomodel->FetchDataToInput(elements,LoadingforceXEnum);
    213219        iomodel->FetchDataToInput(elements,LoadingforceYEnum);
    214         iomodel->FetchDataToInput(elements,DamageDEnum);
     220        if(isdamage)iomodel->FetchDataToInput(elements,DamageDEnum);
    215221
    216222        if(iomodel->domaintype==Domain3DEnum){

  • issm/trunk-jpl/src/c/classes/Materials/Matice.cpp

    r17619 r17757  
    4848        this->element=NULL;
    4949
     50         /*Other perporties*/
     51   int    materialtype;
     52   iomodel->Constant(&materialtype,MaterialsEnum);
     53   if(materialtype==MatdamageiceEnum) this->isdamage = true;
     54   else if(materialtype==MaticeEnum) this->isdamage = false;
     55   else _error_("Material type not recognized");
     56
    5057        return;
    5158
    … …  
    199206IssmDouble Matice::GetD(){
    200207
     208        _assert_(this->isdamage);
    201209        /*Output*/
    202210        IssmDouble D;
    203 
    204211        element->inputs->GetInputAverage(&D,DamageDEnum);
    205212        return D;
    … …  
    209216IssmDouble Matice::GetDbar(){
    210217
     218        _assert_(this->isdamage);
    211219        /*Output*/
    212220        IssmDouble Dbar;
    … …  
    233241
    234242        /*Intermediary: */
    235         IssmDouble B,D,n;
     243        IssmDouble B,D=0.,n;
    236244
    237245        /*Get B and n*/
    238246        B=GetB(); _assert_(B>0.);
    239247        n=GetN(); _assert_(n>0.);
    240         D=GetD(); _assert_(D>=0. && D<1.);
     248        if(this->isdamage){
     249                D=GetD();
     250                _assert_(D>=0. && D<1.);
     251        }
    241252
    242253        if (n==1.){
    … …  
    283294
    284295        /*Intermediary: */
    285         IssmDouble B,D,n;
     296        IssmDouble B,D=0.,n;
    286297
    287298        /*Get B and n*/
    288299        B=GetBbar(); _assert_(B>0.);
    289300        n=GetN();    _assert_(n>0.);
    290         D=GetDbar(); _assert_(D>=0. && D<1.);
     301        if(this->isdamage){
     302                D=GetDbar();
     303                _assert_(D>=0. && D<1.);
     304        }
    291305
    292306        if (n==1.){
    … …  
    334348        /*Intermediary value A and exponent e: */
    335349        IssmDouble A,e;
    336         IssmDouble D,n;
     350        IssmDouble D=0.,n;
    337351
    338352        /*Get D and n*/
    339         D=GetDbar();
     353        if(this->isdamage){
     354                D=GetDbar(); /* GetD()? */
     355                _assert_(D>=0. && D<1.);
     356        }
    340357        n=GetN();
    341358

  • issm/trunk-jpl/src/c/classes/Materials/Matice.h

    r17281 r17757  
    2424        private:
    2525                int      mid;
     26                bool     isdamage;
    2627                Hook    *helement;
    2728                Element *element;

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

    r17692 r17757  
    5252                        iomodel->FetchDataToInput(elements,MaterialsRheologyBEnum);
    5353                        iomodel->FetchDataToInput(elements,MaterialsRheologyNEnum);
     54                        for (i=0;i<iomodel->numberofelements;i++) if(iomodel->my_elements[i]) materials->AddObject(new Matice(i+1,i,iomodel));
     55                        switch(iomodel->domaindim){
     56                                case 2:
     57                                        elements->InputDuplicate(MaterialsRheologyBEnum,MaterialsRheologyBbarEnum);
     58                                        if(dakota_analysis) elements->InputDuplicate(MaterialsRheologyBbarEnum,QmuMaterialsRheologyBEnum);
     59                                        break;
     60                                case 3:
     61                                        if(dakota_analysis) elements->InputDuplicate(MaterialsRheologyBEnum,QmuMaterialsRheologyBEnum);
     62                                        break;
     63                                default:
     64                                        _error_("Mesh not supported yet");
     65                        }
     66                        break;
     67                case MatdamageiceEnum:
     68                        iomodel->FetchDataToInput(elements,MaterialsRheologyBEnum);
     69                        iomodel->FetchDataToInput(elements,MaterialsRheologyNEnum);
    5470                        iomodel->FetchDataToInput(elements,DamageDEnum);
    5571                        for (i=0;i<iomodel->numberofelements;i++) if(iomodel->my_elements[i]) materials->AddObject(new Matice(i+1,i,iomodel));

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

    r17734 r17757  
    1919
    2020        int         i,j,m,k;
    21         int         numoutputs,domaintype,smb_model;
     21        int         numoutputs,domaintype,materialtype,smb_model;
    2222        char**      requestedoutputs = NULL;
    2323        IssmDouble  time;
    … …  
    149149        iomodel->DeleteData(&requestedoutputs,numoutputs,SteadystateRequestedOutputsEnum);
    150150
    151         iomodel->FetchData(&requestedoutputs,&numoutputs,DamageEvolutionRequestedOutputsEnum);
    152         parameters->AddObject(new IntParam(DamageEvolutionNumRequestedOutputsEnum,numoutputs));
    153         if(numoutputs)parameters->AddObject(new StringArrayParam(DamageEvolutionRequestedOutputsEnum,requestedoutputs,numoutputs));
    154         iomodel->DeleteData(&requestedoutputs,numoutputs,DamageEvolutionRequestedOutputsEnum);
     151        iomodel->Constant(&materialtype,MaterialsEnum);
     152        if(materialtype==MatdamageiceEnum){
     153                iomodel->FetchData(&requestedoutputs,&numoutputs,DamageEvolutionRequestedOutputsEnum);
     154                parameters->AddObject(new IntParam(DamageEvolutionNumRequestedOutputsEnum,numoutputs));
     155                if(numoutputs)parameters->AddObject(new StringArrayParam(DamageEvolutionRequestedOutputsEnum,requestedoutputs,numoutputs));
     156                iomodel->DeleteData(&requestedoutputs,numoutputs,DamageEvolutionRequestedOutputsEnum);
     157        }
    155158
    156159        /*Deal with mass flux segments: {{{*/

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

    r17737 r17757  
    8282                if(solution_enum==TransientSolutionEnum && analysis_enum==HydrologyDCEfficientAnalysisEnum && ishydrology==false) continue;
    8383                if(solution_enum==TransientSolutionEnum && analysis_enum==HydrologyDCInefficientAnalysisEnum && ishydrology==false) continue;
     84                if(solution_enum==TransientSolutionEnum && analysis_enum==DamageEvolutionAnalysisEnum && isdamage==false) continue;
    8485
    8586                if(VerboseMProcessor()) _printf0_("   creating datasets for analysis " << EnumToStringx(analysis_enum) << "\n");

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

    r17742 r17757  
    173173        DamageC3Enum,
    174174        DamageC4Enum,
     175        DamageElementinterpEnum,
    175176        DamageHealingEnum,
    176177        DamageStressThresholdEnum,
    … …  
    416417        TransientParamEnum,
    417418        MaticeEnum,
     419        MatdamageiceEnum,
    418420        MatparEnum,
    419421        NodeEnum,

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

    r17742 r17757  
    181181                case DamageC3Enum : return "DamageC3";
    182182                case DamageC4Enum : return "DamageC4";
     183                case DamageElementinterpEnum : return "DamageElementinterp";
    183184                case DamageHealingEnum : return "DamageHealing";
    184185                case DamageStressThresholdEnum : return "DamageStressThreshold";
    … …  
    411412                case TransientParamEnum : return "TransientParam";
    412413                case MaticeEnum : return "Matice";
     414                case MatdamageiceEnum : return "Matdamageice";
    413415                case MatparEnum : return "Matpar";
    414416                case NodeEnum : return "Node";

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

    r17742 r17757  
    184184              else if (strcmp(name,"DamageC3")==0) return DamageC3Enum;
    185185              else if (strcmp(name,"DamageC4")==0) return DamageC4Enum;
     186              else if (strcmp(name,"DamageElementinterp")==0) return DamageElementinterpEnum;
    186187              else if (strcmp(name,"DamageHealing")==0) return DamageHealingEnum;
    187188              else if (strcmp(name,"DamageStressThreshold")==0) return DamageStressThresholdEnum;
    … …  
    259260              else if (strcmp(name,"SettingsWaitonlock")==0) return SettingsWaitonlockEnum;
    260261              else if (strcmp(name,"DebugProfiling")==0) return DebugProfilingEnum;
    261               else if (strcmp(name,"ProfilingCurrentMem")==0) return ProfilingCurrentMemEnum;
    262262         else stage=3;
    263263   }
    264264   if(stage==3){
    265               if (strcmp(name,"ProfilingCurrentFlops")==0) return ProfilingCurrentFlopsEnum;
     265              if (strcmp(name,"ProfilingCurrentMem")==0) return ProfilingCurrentMemEnum;
     266              else if (strcmp(name,"ProfilingCurrentFlops")==0) return ProfilingCurrentFlopsEnum;
    266267              else if (strcmp(name,"ProfilingSolutionTime")==0) return ProfilingSolutionTimeEnum;
    267268              else if (strcmp(name,"MaxIterationConvergenceFlag")==0) return MaxIterationConvergenceFlagEnum;
    … …  
    382383              else if (strcmp(name,"HOFSApproximation")==0) return HOFSApproximationEnum;
    383384              else if (strcmp(name,"FSApproximation")==0) return FSApproximationEnum;
    384               else if (strcmp(name,"FSvelocity")==0) return FSvelocityEnum;
    385385         else stage=4;
    386386   }
    387387   if(stage==4){
    388               if (strcmp(name,"FSpressure")==0) return FSpressureEnum;
     388              if (strcmp(name,"FSvelocity")==0) return FSvelocityEnum;
     389              else if (strcmp(name,"FSpressure")==0) return FSpressureEnum;
    389390              else if (strcmp(name,"Constraints")==0) return ConstraintsEnum;
    390391              else if (strcmp(name,"Loads")==0) return LoadsEnum;
    … …  
    420421              else if (strcmp(name,"TransientParam")==0) return TransientParamEnum;
    421422              else if (strcmp(name,"Matice")==0) return MaticeEnum;
     423              else if (strcmp(name,"Matdamageice")==0) return MatdamageiceEnum;
    422424              else if (strcmp(name,"Matpar")==0) return MatparEnum;
    423425              else if (strcmp(name,"Node")==0) return NodeEnum;
    … …  
    504506              else if (strcmp(name,"Vel")==0) return VelEnum;
    505507              else if (strcmp(name,"Velocity")==0) return VelocityEnum;
    506               else if (strcmp(name,"VxAverage")==0) return VxAverageEnum;
    507               else if (strcmp(name,"Vx")==0) return VxEnum;
    508508         else stage=5;
    509509   }
    510510   if(stage==5){
    511               if (strcmp(name,"VxPicard")==0) return VxPicardEnum;
     511              if (strcmp(name,"VxAverage")==0) return VxAverageEnum;
     512              else if (strcmp(name,"Vx")==0) return VxEnum;
     513              else if (strcmp(name,"VxPicard")==0) return VxPicardEnum;
    512514              else if (strcmp(name,"VyAverage")==0) return VyAverageEnum;
    513515              else if (strcmp(name,"Vy")==0) return VyEnum;
    … …  
    627629              else if (strcmp(name,"SubelementMigration2")==0) return SubelementMigration2Enum;
    628630              else if (strcmp(name,"Contact")==0) return ContactEnum;
    629               else if (strcmp(name,"MaskGroundediceLevelset")==0) return MaskGroundediceLevelsetEnum;
    630               else if (strcmp(name,"QmuMaskGroundediceLevelset")==0) return QmuMaskGroundediceLevelsetEnum;
    631631         else stage=6;
    632632   }
    633633   if(stage==6){
    634               if (strcmp(name,"GaussSeg")==0) return GaussSegEnum;
     634              if (strcmp(name,"MaskGroundediceLevelset")==0) return MaskGroundediceLevelsetEnum;
     635              else if (strcmp(name,"QmuMaskGroundediceLevelset")==0) return QmuMaskGroundediceLevelsetEnum;
     636              else if (strcmp(name,"GaussSeg")==0) return GaussSegEnum;
    635637              else if (strcmp(name,"GaussTria")==0) return GaussTriaEnum;
    636638              else if (strcmp(name,"GaussTetra")==0) return GaussTetraEnum;

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

    r17745 r17757  
    1414                %numerical
    1515                stabilization       = NaN;
     16                maxiter             = NaN;
     17                elementinterp       = '';
    1618                penalty_threshold   = NaN;
    17                 maxiter             = NaN;
    1819                penalty_lock        = NaN;
    1920                penalty_factor      = NaN;
    … …  
    110111                        obj.maxiter=100;
    111112
     113                        %finite element interpolation
     114                        obj.elementinterp='P1';
     115
    112116                        %factor used to compute the values of the penalties: kappa=max(stiffness matrix)*10^penalty_factor
    113117                        obj.penalty_factor=3;
    … …  
    134138                function md = checkconsistency(obj,md,solution,analyses) % {{{
    135139                       
    136                         md = checkfield(md,'fieldname','damage.D','>=',0,'<=',obj.max_damage,'size',[md.mesh.numberofvertices 1]);
    137                         md = checkfield(md,'fieldname','damage.max_damage','<',1,'>=',0);
    138140                        md = checkfield(md,'fieldname','damage.law','values',{'undamaged','pralong'});
    139                         md = checkfield(md,'fieldname','damage.spcdamage','forcing',1);
    140                        
    141                         md = checkfield(md,'fieldname','damage.stabilization','numel',[1],'values',[0 1 2]);
    142                         md = checkfield(md,'fieldname','damage.maxiter','>=0',0);
    143                         md = checkfield(md,'fieldname','damage.penalty_factor','>=0',0);
    144                         md = checkfield(md,'fieldname','damage.penalty_lock','>=0',0);
    145                         md = checkfield(md,'fieldname','damage.penalty_threshold','>=0',0);
     141                        if ~strcmpi(obj.law,'undamaged'),
     142                                md = checkfield(md,'fieldname','damage.D','>=',0,'<=',obj.max_damage,'size',[md.mesh.numberofvertices 1]);
     143                                md = checkfield(md,'fieldname','damage.max_damage','<',1,'>=',0);
     144                                md = checkfield(md,'fieldname','damage.spcdamage','forcing',1);
     145                               
     146                                md = checkfield(md,'fieldname','damage.stabilization','numel',[1],'values',[0 1 2]);
     147                                md = checkfield(md,'fieldname','damage.maxiter','>=0',0);
     148                                md = checkfield(md,'fieldname','damage.elementinterp','values',{'P1','P2'});
     149                                md = checkfield(md,'fieldname','damage.penalty_factor','>=0',0);
     150                                md = checkfield(md,'fieldname','damage.penalty_lock','>=0',0);
     151                                md = checkfield(md,'fieldname','damage.penalty_threshold','>=0',0);
     152                        end
    146153
    147154                        if strcmpi(obj.law,'pralong'),
    … …  
    169176         else
    170177            list = {'DamageD'};
    171          end
     178                        end
    172179                end % }}}
    173180                function disp(obj) % {{{
    174181                        disp(sprintf('   Damage:\n'));
    175182
    176                         fielddisplay(obj,'D','damage tensor (scalar)');
    177183                        fielddisplay(obj,'law','damage law (string) from {''undamaged'',''pralong''}');
    178                         fielddisplay(obj,'spcdamage','damage constraints (NaN means no constraint)');
    179                         fielddisplay(obj,'max_damage','maximum possible damage (0<=max_damage<1)');
    180                        
    181                         fielddisplay(obj,'stabilization','0: no, 1: artificial_diffusivity, 2: SUPG');
    182                         fielddisplay(obj,'maxiter','maximum number of non linear iterations');
    183                         fielddisplay(obj,'penalty_lock','stabilize unstable damage constraints that keep zigzagging after n iteration (default is 0, no stabilization)');
    184                         fielddisplay(obj,'penalty_threshold','threshold to declare convergence of damage evolution solution (default is 0)');
    185                         fielddisplay(obj,'penalty_factor','scaling exponent (default is 3)');
     184                        if ~strcmpi(obj.law,'undamaged'),
     185                                fielddisplay(obj,'D','damage tensor (scalar)');
     186                                fielddisplay(obj,'spcdamage','damage constraints (NaN means no constraint)');
     187                                fielddisplay(obj,'max_damage','maximum possible damage (0<=max_damage<1)');
     188                               
     189                                fielddisplay(obj,'stabilization','0: no, 1: artificial_diffusivity, 2: SUPG');
     190                                fielddisplay(obj,'maxiter','maximum number of non linear iterations');
     191                                fielddisplay(obj,'elementinterp','interpolation scheme for finite elements {''P1'',''P2''}');
     192                                fielddisplay(obj,'penalty_lock','stabilize unstable damage constraints that keep zigzagging after n iteration (default is 0, no stabilization)');
     193                                fielddisplay(obj,'penalty_threshold','threshold to declare convergence of damage evolution solution (default is 0)');
     194                                fielddisplay(obj,'penalty_factor','scaling exponent (default is 3)');
     195                        end
    186196
    187197                        if strcmpi(obj.law,'pralong'),
    … …  
    199209                function marshall(obj,md,fid) % {{{
    200210               
    201                         WriteData(fid,'object',obj,'fieldname','D','format','DoubleMat','mattype',1);
    202211                        WriteData(fid,'object',obj,'fieldname','law','format','String');
    203                         WriteData(fid,'object',obj,'fieldname','spcdamage','format','DoubleMat','mattype',1,'forcinglength',md.mesh.numberofvertices+1);
    204                         WriteData(fid,'object',obj,'fieldname','max_damage','format','Double');
    205 
    206                         WriteData(fid,'object',obj,'fieldname','stabilization','format','Integer');
    207                         WriteData(fid,'object',obj,'fieldname','penalty_threshold','format','Integer');
    208                         WriteData(fid,'object',obj,'fieldname','maxiter','format','Integer');
    209                         WriteData(fid,'object',obj,'fieldname','penalty_lock','format','Integer');
    210                         WriteData(fid,'object',obj,'fieldname','penalty_factor','format','Double');
     212                        if ~strcmpi(obj.law,'undamaged'),
     213                                WriteData(fid,'object',obj,'fieldname','D','format','DoubleMat','mattype',1);
     214                                WriteData(fid,'object',obj,'fieldname','spcdamage','format','DoubleMat','mattype',1,'forcinglength',md.mesh.numberofvertices+1);
     215                                WriteData(fid,'object',obj,'fieldname','max_damage','format','Double');
     216
     217                                WriteData(fid,'object',obj,'fieldname','stabilization','format','Integer');
     218                                WriteData(fid,'object',obj,'fieldname','maxiter','format','Integer');
     219                                WriteData(fid,'object',obj,'fieldname','elementinterp','format','String');
     220                                WriteData(fid,'object',obj,'fieldname','penalty_threshold','format','Integer');
     221                                WriteData(fid,'object',obj,'fieldname','penalty_lock','format','Integer');
     222                                WriteData(fid,'object',obj,'fieldname','penalty_factor','format','Double');
     223                        end
    211224       
    212225                        if strcmpi(obj.law,'pralong'),
    … …  
    227240                                outputs      = [outputs defaultoutputs(obj,md)]; %add defaults
    228241                        end
    229                         WriteData(fid,'data',outputs,'enum',DamageEvolutionRequestedOutputsEnum,'format','StringArray');
     242                        if ~strcmpi(obj.law,'undamaged'),
     243                                WriteData(fid,'data',outputs,'enum',DamageEvolutionRequestedOutputsEnum,'format','StringArray');
     244                        end
    230245
    231246                end % }}}

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

    r17754 r17757  
    3030
    3131                balancethickness = 0;
    32                 stressbalance       = 0;
     32                stressbalance    = 0;
    3333                groundingline    = 0;
    3434                hydrology        = 0;
    35                 masstransport       = 0;
     35                masstransport    = 0;
    3636                thermal          = 0;
    3737                steadystate      = 0;

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

    r17742 r17757  
    173173def DamageC3Enum(): return StringToEnum("DamageC3")[0]
    174174def DamageC4Enum(): return StringToEnum("DamageC4")[0]
     175def DamageElementinterpEnum(): return StringToEnum("DamageElementinterp")[0]
    175176def DamageHealingEnum(): return StringToEnum("DamageHealing")[0]
    176177def DamageStressThresholdEnum(): return StringToEnum("DamageStressThreshold")[0]
    … …  
    403404def TransientParamEnum(): return StringToEnum("TransientParam")[0]
    404405def MaticeEnum(): return StringToEnum("Matice")[0]
     406def MatdamageiceEnum(): return StringToEnum("Matdamageice")[0]
    405407def MatparEnum(): return StringToEnum("Matpar")[0]
    406408def NodeEnum(): return StringToEnum("Node")[0]

  • issm/trunk-jpl/src/m/mech/mechanicalproperties.m

    r17686 r17757  
    116116stress.principalvalue2=valuesstress(:,2);
    117117stress.principalaxis2=directionsstress(:,3:4);
    118 stress.effectivevalue=1/sqrt(2)*sqrt(stress.xx.^2+stress.yy.^2+2*stress.xy.^2); % effective shear stress
     118stress.effectivevalue=1/sqrt(2)*sqrt(stress.xx.^2+stress.yy.^2+2*stress.xy.^2); % effective shear stress (sqrt(J2))
    119119md.results.stress=stress;
    120120
Note: See TracChangeset for help on using the changeset viewer.