Index: /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp (revision 17756) +++ /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp (revision 17757) @@ -131,4 +131,5 @@ bool dakota_analysis; bool islevelset; + bool isdamage; /*Fetch constants needed: */ @@ -154,4 +155,9 @@ else iscoupling = false; + + /*is damage mechanics being used?*/ + if(materials_type==MaticeEnum) isdamage = false; + else if(materials_type==MatdamageiceEnum) isdamage = true; + else _error_("Material type not recognized"); /*Get finite element type*/ @@ -212,5 +218,5 @@ iomodel->FetchDataToInput(elements,LoadingforceXEnum); iomodel->FetchDataToInput(elements,LoadingforceYEnum); - iomodel->FetchDataToInput(elements,DamageDEnum); + if(isdamage)iomodel->FetchDataToInput(elements,DamageDEnum); if(iomodel->domaintype==Domain3DEnum){ Index: /issm/trunk-jpl/src/c/classes/Materials/Matice.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Materials/Matice.cpp (revision 17756) +++ /issm/trunk-jpl/src/c/classes/Materials/Matice.cpp (revision 17757) @@ -48,4 +48,11 @@ this->element=NULL; + /*Other perporties*/ + int materialtype; + iomodel->Constant(&materialtype,MaterialsEnum); + if(materialtype==MatdamageiceEnum) this->isdamage = true; + else if(materialtype==MaticeEnum) this->isdamage = false; + else _error_("Material type not recognized"); + return; @@ -199,7 +206,7 @@ IssmDouble Matice::GetD(){ + _assert_(this->isdamage); /*Output*/ IssmDouble D; - element->inputs->GetInputAverage(&D,DamageDEnum); return D; @@ -209,4 +216,5 @@ IssmDouble Matice::GetDbar(){ + _assert_(this->isdamage); /*Output*/ IssmDouble Dbar; @@ -233,10 +241,13 @@ /*Intermediary: */ - IssmDouble B,D,n; + IssmDouble B,D=0.,n; /*Get B and n*/ B=GetB(); _assert_(B>0.); n=GetN(); _assert_(n>0.); - D=GetD(); _assert_(D>=0. && D<1.); + if(this->isdamage){ + D=GetD(); + _assert_(D>=0. && D<1.); + } if (n==1.){ @@ -283,10 +294,13 @@ /*Intermediary: */ - IssmDouble B,D,n; + IssmDouble B,D=0.,n; /*Get B and n*/ B=GetBbar(); _assert_(B>0.); n=GetN(); _assert_(n>0.); - D=GetDbar(); _assert_(D>=0. && D<1.); + if(this->isdamage){ + D=GetDbar(); + _assert_(D>=0. && D<1.); + } if (n==1.){ @@ -334,8 +348,11 @@ /*Intermediary value A and exponent e: */ IssmDouble A,e; - IssmDouble D,n; + IssmDouble D=0.,n; /*Get D and n*/ - D=GetDbar(); + if(this->isdamage){ + D=GetDbar(); /* GetD()? */ + _assert_(D>=0. && D<1.); + } n=GetN(); Index: /issm/trunk-jpl/src/c/classes/Materials/Matice.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Materials/Matice.h (revision 17756) +++ /issm/trunk-jpl/src/c/classes/Materials/Matice.h (revision 17757) @@ -24,4 +24,5 @@ private: int mid; + bool isdamage; Hook *helement; Element *element; Index: /issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateElementsVerticesAndMaterials.cpp =================================================================== --- /issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateElementsVerticesAndMaterials.cpp (revision 17756) +++ /issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateElementsVerticesAndMaterials.cpp (revision 17757) @@ -52,4 +52,20 @@ iomodel->FetchDataToInput(elements,MaterialsRheologyBEnum); iomodel->FetchDataToInput(elements,MaterialsRheologyNEnum); + for (i=0;inumberofelements;i++) if(iomodel->my_elements[i]) materials->AddObject(new Matice(i+1,i,iomodel)); + switch(iomodel->domaindim){ + case 2: + elements->InputDuplicate(MaterialsRheologyBEnum,MaterialsRheologyBbarEnum); + if(dakota_analysis) elements->InputDuplicate(MaterialsRheologyBbarEnum,QmuMaterialsRheologyBEnum); + break; + case 3: + if(dakota_analysis) elements->InputDuplicate(MaterialsRheologyBEnum,QmuMaterialsRheologyBEnum); + break; + default: + _error_("Mesh not supported yet"); + } + break; + case MatdamageiceEnum: + iomodel->FetchDataToInput(elements,MaterialsRheologyBEnum); + iomodel->FetchDataToInput(elements,MaterialsRheologyNEnum); iomodel->FetchDataToInput(elements,DamageDEnum); for (i=0;inumberofelements;i++) if(iomodel->my_elements[i]) materials->AddObject(new Matice(i+1,i,iomodel)); Index: /issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateParameters.cpp =================================================================== --- /issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateParameters.cpp (revision 17756) +++ /issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateParameters.cpp (revision 17757) @@ -19,5 +19,5 @@ int i,j,m,k; - int numoutputs,domaintype,smb_model; + int numoutputs,domaintype,materialtype,smb_model; char** requestedoutputs = NULL; IssmDouble time; @@ -149,8 +149,11 @@ iomodel->DeleteData(&requestedoutputs,numoutputs,SteadystateRequestedOutputsEnum); - iomodel->FetchData(&requestedoutputs,&numoutputs,DamageEvolutionRequestedOutputsEnum); - parameters->AddObject(new IntParam(DamageEvolutionNumRequestedOutputsEnum,numoutputs)); - if(numoutputs)parameters->AddObject(new StringArrayParam(DamageEvolutionRequestedOutputsEnum,requestedoutputs,numoutputs)); - iomodel->DeleteData(&requestedoutputs,numoutputs,DamageEvolutionRequestedOutputsEnum); + iomodel->Constant(&materialtype,MaterialsEnum); + if(materialtype==MatdamageiceEnum){ + iomodel->FetchData(&requestedoutputs,&numoutputs,DamageEvolutionRequestedOutputsEnum); + parameters->AddObject(new IntParam(DamageEvolutionNumRequestedOutputsEnum,numoutputs)); + if(numoutputs)parameters->AddObject(new StringArrayParam(DamageEvolutionRequestedOutputsEnum,requestedoutputs,numoutputs)); + iomodel->DeleteData(&requestedoutputs,numoutputs,DamageEvolutionRequestedOutputsEnum); + } /*Deal with mass flux segments: {{{*/ Index: /issm/trunk-jpl/src/c/modules/ModelProcessorx/ModelProcessorx.cpp =================================================================== --- /issm/trunk-jpl/src/c/modules/ModelProcessorx/ModelProcessorx.cpp (revision 17756) +++ /issm/trunk-jpl/src/c/modules/ModelProcessorx/ModelProcessorx.cpp (revision 17757) @@ -82,4 +82,5 @@ if(solution_enum==TransientSolutionEnum && analysis_enum==HydrologyDCEfficientAnalysisEnum && ishydrology==false) continue; if(solution_enum==TransientSolutionEnum && analysis_enum==HydrologyDCInefficientAnalysisEnum && ishydrology==false) continue; + if(solution_enum==TransientSolutionEnum && analysis_enum==DamageEvolutionAnalysisEnum && isdamage==false) continue; if(VerboseMProcessor()) _printf0_(" creating datasets for analysis " << EnumToStringx(analysis_enum) << "\n"); Index: /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h (revision 17756) +++ /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h (revision 17757) @@ -173,4 +173,5 @@ DamageC3Enum, DamageC4Enum, + DamageElementinterpEnum, DamageHealingEnum, DamageStressThresholdEnum, @@ -416,4 +417,5 @@ TransientParamEnum, MaticeEnum, + MatdamageiceEnum, MatparEnum, NodeEnum, Index: /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp (revision 17756) +++ /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp (revision 17757) @@ -181,4 +181,5 @@ case DamageC3Enum : return "DamageC3"; case DamageC4Enum : return "DamageC4"; + case DamageElementinterpEnum : return "DamageElementinterp"; case DamageHealingEnum : return "DamageHealing"; case DamageStressThresholdEnum : return "DamageStressThreshold"; @@ -411,4 +412,5 @@ case TransientParamEnum : return "TransientParam"; case MaticeEnum : return "Matice"; + case MatdamageiceEnum : return "Matdamageice"; case MatparEnum : return "Matpar"; case NodeEnum : return "Node"; Index: /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp (revision 17756) +++ /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp (revision 17757) @@ -184,4 +184,5 @@ else if (strcmp(name,"DamageC3")==0) return DamageC3Enum; else if (strcmp(name,"DamageC4")==0) return DamageC4Enum; + else if (strcmp(name,"DamageElementinterp")==0) return DamageElementinterpEnum; else if (strcmp(name,"DamageHealing")==0) return DamageHealingEnum; else if (strcmp(name,"DamageStressThreshold")==0) return DamageStressThresholdEnum; @@ -259,9 +260,9 @@ else if (strcmp(name,"SettingsWaitonlock")==0) return SettingsWaitonlockEnum; else if (strcmp(name,"DebugProfiling")==0) return DebugProfilingEnum; - else if (strcmp(name,"ProfilingCurrentMem")==0) return ProfilingCurrentMemEnum; else stage=3; } if(stage==3){ - if (strcmp(name,"ProfilingCurrentFlops")==0) return ProfilingCurrentFlopsEnum; + if (strcmp(name,"ProfilingCurrentMem")==0) return ProfilingCurrentMemEnum; + else if (strcmp(name,"ProfilingCurrentFlops")==0) return ProfilingCurrentFlopsEnum; else if (strcmp(name,"ProfilingSolutionTime")==0) return ProfilingSolutionTimeEnum; else if (strcmp(name,"MaxIterationConvergenceFlag")==0) return MaxIterationConvergenceFlagEnum; @@ -382,9 +383,9 @@ else if (strcmp(name,"HOFSApproximation")==0) return HOFSApproximationEnum; else if (strcmp(name,"FSApproximation")==0) return FSApproximationEnum; - else if (strcmp(name,"FSvelocity")==0) return FSvelocityEnum; else stage=4; } if(stage==4){ - if (strcmp(name,"FSpressure")==0) return FSpressureEnum; + if (strcmp(name,"FSvelocity")==0) return FSvelocityEnum; + else if (strcmp(name,"FSpressure")==0) return FSpressureEnum; else if (strcmp(name,"Constraints")==0) return ConstraintsEnum; else if (strcmp(name,"Loads")==0) return LoadsEnum; @@ -420,4 +421,5 @@ else if (strcmp(name,"TransientParam")==0) return TransientParamEnum; else if (strcmp(name,"Matice")==0) return MaticeEnum; + else if (strcmp(name,"Matdamageice")==0) return MatdamageiceEnum; else if (strcmp(name,"Matpar")==0) return MatparEnum; else if (strcmp(name,"Node")==0) return NodeEnum; @@ -504,10 +506,10 @@ else if (strcmp(name,"Vel")==0) return VelEnum; else if (strcmp(name,"Velocity")==0) return VelocityEnum; - else if (strcmp(name,"VxAverage")==0) return VxAverageEnum; - else if (strcmp(name,"Vx")==0) return VxEnum; else stage=5; } if(stage==5){ - if (strcmp(name,"VxPicard")==0) return VxPicardEnum; + if (strcmp(name,"VxAverage")==0) return VxAverageEnum; + else if (strcmp(name,"Vx")==0) return VxEnum; + else if (strcmp(name,"VxPicard")==0) return VxPicardEnum; else if (strcmp(name,"VyAverage")==0) return VyAverageEnum; else if (strcmp(name,"Vy")==0) return VyEnum; @@ -627,10 +629,10 @@ else if (strcmp(name,"SubelementMigration2")==0) return SubelementMigration2Enum; else if (strcmp(name,"Contact")==0) return ContactEnum; - else if (strcmp(name,"MaskGroundediceLevelset")==0) return MaskGroundediceLevelsetEnum; - else if (strcmp(name,"QmuMaskGroundediceLevelset")==0) return QmuMaskGroundediceLevelsetEnum; else stage=6; } if(stage==6){ - if (strcmp(name,"GaussSeg")==0) return GaussSegEnum; + if (strcmp(name,"MaskGroundediceLevelset")==0) return MaskGroundediceLevelsetEnum; + else if (strcmp(name,"QmuMaskGroundediceLevelset")==0) return QmuMaskGroundediceLevelsetEnum; + else if (strcmp(name,"GaussSeg")==0) return GaussSegEnum; else if (strcmp(name,"GaussTria")==0) return GaussTriaEnum; else if (strcmp(name,"GaussTetra")==0) return GaussTetraEnum; Index: /issm/trunk-jpl/src/m/classes/damage.m =================================================================== --- /issm/trunk-jpl/src/m/classes/damage.m (revision 17756) +++ /issm/trunk-jpl/src/m/classes/damage.m (revision 17757) @@ -14,6 +14,7 @@ %numerical stabilization = NaN; + maxiter = NaN; + elementinterp = ''; penalty_threshold = NaN; - maxiter = NaN; penalty_lock = NaN; penalty_factor = NaN; @@ -110,4 +111,7 @@ obj.maxiter=100; + %finite element interpolation + obj.elementinterp='P1'; + %factor used to compute the values of the penalties: kappa=max(stiffness matrix)*10^penalty_factor obj.penalty_factor=3; @@ -134,14 +138,17 @@ function md = checkconsistency(obj,md,solution,analyses) % {{{ - md = checkfield(md,'fieldname','damage.D','>=',0,'<=',obj.max_damage,'size',[md.mesh.numberofvertices 1]); - md = checkfield(md,'fieldname','damage.max_damage','<',1,'>=',0); md = checkfield(md,'fieldname','damage.law','values',{'undamaged','pralong'}); - md = checkfield(md,'fieldname','damage.spcdamage','forcing',1); - - md = checkfield(md,'fieldname','damage.stabilization','numel',[1],'values',[0 1 2]); - md = checkfield(md,'fieldname','damage.maxiter','>=0',0); - md = checkfield(md,'fieldname','damage.penalty_factor','>=0',0); - md = checkfield(md,'fieldname','damage.penalty_lock','>=0',0); - md = checkfield(md,'fieldname','damage.penalty_threshold','>=0',0); + if ~strcmpi(obj.law,'undamaged'), + md = checkfield(md,'fieldname','damage.D','>=',0,'<=',obj.max_damage,'size',[md.mesh.numberofvertices 1]); + md = checkfield(md,'fieldname','damage.max_damage','<',1,'>=',0); + md = checkfield(md,'fieldname','damage.spcdamage','forcing',1); + + md = checkfield(md,'fieldname','damage.stabilization','numel',[1],'values',[0 1 2]); + md = checkfield(md,'fieldname','damage.maxiter','>=0',0); + md = checkfield(md,'fieldname','damage.elementinterp','values',{'P1','P2'}); + md = checkfield(md,'fieldname','damage.penalty_factor','>=0',0); + md = checkfield(md,'fieldname','damage.penalty_lock','>=0',0); + md = checkfield(md,'fieldname','damage.penalty_threshold','>=0',0); + end if strcmpi(obj.law,'pralong'), @@ -169,19 +176,22 @@ else list = {'DamageD'}; - end + end end % }}} function disp(obj) % {{{ disp(sprintf(' Damage:\n')); - fielddisplay(obj,'D','damage tensor (scalar)'); fielddisplay(obj,'law','damage law (string) from {''undamaged'',''pralong''}'); - fielddisplay(obj,'spcdamage','damage constraints (NaN means no constraint)'); - fielddisplay(obj,'max_damage','maximum possible damage (0<=max_damage<1)'); - - fielddisplay(obj,'stabilization','0: no, 1: artificial_diffusivity, 2: SUPG'); - fielddisplay(obj,'maxiter','maximum number of non linear iterations'); - fielddisplay(obj,'penalty_lock','stabilize unstable damage constraints that keep zigzagging after n iteration (default is 0, no stabilization)'); - fielddisplay(obj,'penalty_threshold','threshold to declare convergence of damage evolution solution (default is 0)'); - fielddisplay(obj,'penalty_factor','scaling exponent (default is 3)'); + if ~strcmpi(obj.law,'undamaged'), + fielddisplay(obj,'D','damage tensor (scalar)'); + fielddisplay(obj,'spcdamage','damage constraints (NaN means no constraint)'); + fielddisplay(obj,'max_damage','maximum possible damage (0<=max_damage<1)'); + + fielddisplay(obj,'stabilization','0: no, 1: artificial_diffusivity, 2: SUPG'); + fielddisplay(obj,'maxiter','maximum number of non linear iterations'); + fielddisplay(obj,'elementinterp','interpolation scheme for finite elements {''P1'',''P2''}'); + fielddisplay(obj,'penalty_lock','stabilize unstable damage constraints that keep zigzagging after n iteration (default is 0, no stabilization)'); + fielddisplay(obj,'penalty_threshold','threshold to declare convergence of damage evolution solution (default is 0)'); + fielddisplay(obj,'penalty_factor','scaling exponent (default is 3)'); + end if strcmpi(obj.law,'pralong'), @@ -199,14 +209,17 @@ function marshall(obj,md,fid) % {{{ - WriteData(fid,'object',obj,'fieldname','D','format','DoubleMat','mattype',1); WriteData(fid,'object',obj,'fieldname','law','format','String'); - WriteData(fid,'object',obj,'fieldname','spcdamage','format','DoubleMat','mattype',1,'forcinglength',md.mesh.numberofvertices+1); - WriteData(fid,'object',obj,'fieldname','max_damage','format','Double'); - - WriteData(fid,'object',obj,'fieldname','stabilization','format','Integer'); - WriteData(fid,'object',obj,'fieldname','penalty_threshold','format','Integer'); - WriteData(fid,'object',obj,'fieldname','maxiter','format','Integer'); - WriteData(fid,'object',obj,'fieldname','penalty_lock','format','Integer'); - WriteData(fid,'object',obj,'fieldname','penalty_factor','format','Double'); + if ~strcmpi(obj.law,'undamaged'), + WriteData(fid,'object',obj,'fieldname','D','format','DoubleMat','mattype',1); + WriteData(fid,'object',obj,'fieldname','spcdamage','format','DoubleMat','mattype',1,'forcinglength',md.mesh.numberofvertices+1); + WriteData(fid,'object',obj,'fieldname','max_damage','format','Double'); + + WriteData(fid,'object',obj,'fieldname','stabilization','format','Integer'); + WriteData(fid,'object',obj,'fieldname','maxiter','format','Integer'); + WriteData(fid,'object',obj,'fieldname','elementinterp','format','String'); + WriteData(fid,'object',obj,'fieldname','penalty_threshold','format','Integer'); + WriteData(fid,'object',obj,'fieldname','penalty_lock','format','Integer'); + WriteData(fid,'object',obj,'fieldname','penalty_factor','format','Double'); + end if strcmpi(obj.law,'pralong'), @@ -227,5 +240,7 @@ outputs = [outputs defaultoutputs(obj,md)]; %add defaults end - WriteData(fid,'data',outputs,'enum',DamageEvolutionRequestedOutputsEnum,'format','StringArray'); + if ~strcmpi(obj.law,'undamaged'), + WriteData(fid,'data',outputs,'enum',DamageEvolutionRequestedOutputsEnum,'format','StringArray'); + end end % }}} Index: /issm/trunk-jpl/src/m/classes/matdamageice.m =================================================================== --- /issm/trunk-jpl/src/m/classes/matdamageice.m (revision 17757) +++ /issm/trunk-jpl/src/m/classes/matdamageice.m (revision 17757) @@ -0,0 +1,197 @@ +%MATICE class definition +% +% Usage: +% matice=matice(); + +classdef matice + properties (SetAccess=public) + rho_ice = 0.; + rho_water = 0.; + rho_freshwater = 0.; + mu_water = 0.; + heatcapacity = 0.; + latentheat = 0.; + thermalconductivity = 0.; + temperateiceconductivity = 0.; + meltingpoint = 0.; + beta = 0.; + mixed_layer_capacity = 0.; + thermal_exchange_velocity = 0.; + rheology_B = NaN; + rheology_n = NaN; + rheology_law = ''; + + %gia: + lithosphere_shear_modulus = 0.; + lithosphere_density = 0.; + mantle_shear_modulus = 0.; + mantle_density = 0.; + + end + methods + function createxml(obj,fid) % {{{ + fprintf(fid, '\n\n'); + fprintf(fid, '\n'); + fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','','
',' ice density [kg/m^3] ',''); + fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','','
',' ocean water density [kg/m^3] ',''); + fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','','
',' fresh water density [kg/m^3] ',''); + + + fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','','
',' water viscosity [N s/m^2] ',''); + fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','','
',' heat capacity [J/kg/K] ',''); + fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','','
',' latent heat of fusion [J/kg] ',''); + + + fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','','
',' ice thermal conductivity [W/m/K] ',''); + fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','','
',' temperate ice thermal conductivity [W/m/K] ',''); + fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','','
',' melting point of ice at 1atm in K ',''); + + + fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','','
',' rate of change of melting point with pressure [K/Pa] ',''); + fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','','
',' mixed layer capacity [W/kg/K] ',''); + fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','','
',' thermal exchange velocity [m/s] ',''); + + + fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','','
',' flow law parameter [Pa/s^(1/n)] ',''); + fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','','
',' Glens flow law exponent ',''); + fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','','
',' law for the temperature dependance of the rheology: "None", "Paterson", "Arrhenius" or "LliboutryDuval" ',''); + + fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','','
',' Lithosphere shear modulus [Pa] ',''); + fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','','
',' Lithosphere density [g/cm^-3] ',''); + fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','','
',' Mantle shear modulus [Pa] ',''); + fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','','
',' Mantle density [g/cm^-3] ',''); + + + end % }}} + function obj = matice(varargin) % {{{ + switch nargin + case 0 + obj=setdefaultparameters(obj); + case 1 + inputstruct=varargin{1}; + list1 = properties('matice'); + list2 = fieldnames(inputstruct); + for i=1:length(list1) + fieldname = list1{i}; + if ismember(fieldname,list2), + obj.(fieldname) = inputstruct.(fieldname); + end + end + otherwise + error('constructor not supported'); + end + end % }}} + function obj = setdefaultparameters(obj) % {{{ + + %ice density (kg/m^3) + obj.rho_ice=917.; + + %ocean water density (kg/m^3) + obj.rho_water=1023.; + + %fresh water density (kg/m^3) + obj.rho_freshwater=1000.; + + %water viscosity (N.s/m^2) + obj.mu_water=0.001787; + + %ice heat capacity cp (J/kg/K) + obj.heatcapacity=2093.; + + %ice latent heat of fusion L (J/kg) + obj.latentheat=3.34*10^5; + + %ice thermal conductivity (W/m/K) + obj.thermalconductivity=2.4; + + %wet ice thermal conductivity (W/m/K) + obj.temperateiceconductivity=.24; + + %the melting point of ice at 1 atmosphere of pressure in K + obj.meltingpoint=273.15; + + %rate of change of melting point with pressure (K/Pa) + obj.beta=9.8*10^-8; + + %mixed layer (ice-water interface) heat capacity (J/kg/K) + obj.mixed_layer_capacity=3974.; + + %thermal exchange velocity (ice-water interface) (m/s) + obj.thermal_exchange_velocity=1.00*10^-4; + + %Rheology law: what is the temperature dependence of B with T + %available: none, paterson and arrhenius + obj.rheology_law='Paterson'; + + % GIA: + obj.lithosphere_shear_modulus = 6.7*10^10; % (Pa) + obj.lithosphere_density = 3.32; % (g/cm^-3) + obj.mantle_shear_modulus = 1.45*10^11; % (Pa) + obj.mantle_density = 3.34; % (g/cm^-3) + + end % }}} + function md = checkconsistency(obj,md,solution,analyses) % {{{ + md = checkfield(md,'fieldname','materials.rho_ice','>',0); + md = checkfield(md,'fieldname','materials.rho_water','>',0); + md = checkfield(md,'fieldname','materials.rho_freshwater','>',0); + md = checkfield(md,'fieldname','materials.mu_water','>',0); + md = checkfield(md,'fieldname','materials.rheology_B','>',0,'size',[md.mesh.numberofvertices 1]); + md = checkfield(md,'fieldname','materials.rheology_n','>',0,'size',[md.mesh.numberofelements 1]); + md = checkfield(md,'fieldname','materials.rheology_law','values',{'None' 'Cuffey' 'Paterson' 'Arrhenius' 'LliboutryDuval'}); + + if ismember(GiaAnalysisEnum(),analyses), + md = checkfield(md,'fieldname','materials.lithosphere_shear_modulus','>',0,'numel',1); + md = checkfield(md,'fieldname','materials.lithosphere_density','>',0,'numel',1); + md = checkfield(md,'fieldname','materials.mantle_shear_modulus','>',0,'numel',1); + md = checkfield(md,'fieldname','materials.mantle_density','>',0,'numel',1); + end + + end % }}} + function disp(obj) % {{{ + disp(sprintf(' Materials:')); + + fielddisplay(obj,'rho_ice','ice density [kg/m^3]'); + fielddisplay(obj,'rho_water','ocean water density [kg/m^3]'); + fielddisplay(obj,'rho_freshwater','fresh water density [kg/m^3]'); + fielddisplay(obj,'mu_water','water viscosity [N s/m^2]'); + fielddisplay(obj,'heatcapacity','heat capacity [J/kg/K]'); + fielddisplay(obj,'thermalconductivity',['ice thermal conductivity [W/m/K]']); + fielddisplay(obj,'temperateiceconductivity','temperate ice thermal conductivity [W/m/K]'); + fielddisplay(obj,'meltingpoint','melting point of ice at 1atm in K'); + fielddisplay(obj,'latentheat','latent heat of fusion [J/kg]'); + fielddisplay(obj,'beta','rate of change of melting point with pressure [K/Pa]'); + fielddisplay(obj,'mixed_layer_capacity','mixed layer capacity [W/kg/K]'); + fielddisplay(obj,'thermal_exchange_velocity','thermal exchange velocity [m/s]'); + fielddisplay(obj,'rheology_B','flow law parameter [Pa/s^(1/n)]'); + fielddisplay(obj,'rheology_n','Glen''s flow law exponent'); + fielddisplay(obj,'rheology_law',['law for the temperature dependance of the rheology: ''None'', ''Cuffey'', ''Paterson'', ''Arrhenius'' or ''LliboutryDuval''']); + fielddisplay(obj,'lithosphere_shear_modulus','Lithosphere shear modulus [Pa]'); + fielddisplay(obj,'lithosphere_density','Lithosphere density [g/cm^-3]'); + fielddisplay(obj,'mantle_shear_modulus','Mantle shear modulus [Pa]'); + fielddisplay(obj,'mantle_density','Mantle density [g/cm^-3]'); + end % }}} + function marshall(obj,md,fid) % {{{ + WriteData(fid,'enum',MaterialsEnum(),'data',MatdamageiceEnum(),'format','Integer'); + WriteData(fid,'object',obj,'class','materials','fieldname','rho_ice','format','Double'); + WriteData(fid,'object',obj,'class','materials','fieldname','rho_water','format','Double'); + WriteData(fid,'object',obj,'class','materials','fieldname','rho_freshwater','format','Double'); + WriteData(fid,'object',obj,'class','materials','fieldname','mu_water','format','Double'); + WriteData(fid,'object',obj,'class','materials','fieldname','heatcapacity','format','Double'); + WriteData(fid,'object',obj,'class','materials','fieldname','latentheat','format','Double'); + WriteData(fid,'object',obj,'class','materials','fieldname','thermalconductivity','format','Double'); + WriteData(fid,'object',obj,'class','materials','fieldname','temperateiceconductivity','format','Double'); + WriteData(fid,'object',obj,'class','materials','fieldname','meltingpoint','format','Double'); + WriteData(fid,'object',obj,'class','materials','fieldname','beta','format','Double'); + WriteData(fid,'object',obj,'class','materials','fieldname','mixed_layer_capacity','format','Double'); + WriteData(fid,'object',obj,'class','materials','fieldname','thermal_exchange_velocity','format','Double'); + WriteData(fid,'object',obj,'class','materials','fieldname','rheology_B','format','DoubleMat','mattype',1); + WriteData(fid,'object',obj,'class','materials','fieldname','rheology_n','format','DoubleMat','mattype',2); + WriteData(fid,'data',StringToEnum(obj.rheology_law),'enum',MaterialsRheologyLawEnum(),'format','Integer'); + + WriteData(fid,'object',obj,'class','materials','fieldname','lithosphere_shear_modulus','format','Double'); + WriteData(fid,'object',obj,'class','materials','fieldname','lithosphere_density','format','Double','scale',10^3); + WriteData(fid,'object',obj,'class','materials','fieldname','mantle_shear_modulus','format','Double'); + WriteData(fid,'object',obj,'class','materials','fieldname','mantle_density','format','Double','scale',10^3); + end % }}} + end +end Index: /issm/trunk-jpl/src/m/classes/model.m =================================================================== --- /issm/trunk-jpl/src/m/classes/model.m (revision 17756) +++ /issm/trunk-jpl/src/m/classes/model.m (revision 17757) @@ -30,8 +30,8 @@ balancethickness = 0; - stressbalance = 0; + stressbalance = 0; groundingline = 0; hydrology = 0; - masstransport = 0; + masstransport = 0; thermal = 0; steadystate = 0; Index: /issm/trunk-jpl/src/m/enum/EnumDefinitions.py =================================================================== --- /issm/trunk-jpl/src/m/enum/EnumDefinitions.py (revision 17756) +++ /issm/trunk-jpl/src/m/enum/EnumDefinitions.py (revision 17757) @@ -173,4 +173,5 @@ def DamageC3Enum(): return StringToEnum("DamageC3")[0] def DamageC4Enum(): return StringToEnum("DamageC4")[0] +def DamageElementinterpEnum(): return StringToEnum("DamageElementinterp")[0] def DamageHealingEnum(): return StringToEnum("DamageHealing")[0] def DamageStressThresholdEnum(): return StringToEnum("DamageStressThreshold")[0] @@ -403,4 +404,5 @@ def TransientParamEnum(): return StringToEnum("TransientParam")[0] def MaticeEnum(): return StringToEnum("Matice")[0] +def MatdamageiceEnum(): return StringToEnum("Matdamageice")[0] def MatparEnum(): return StringToEnum("Matpar")[0] def NodeEnum(): return StringToEnum("Node")[0] Index: /issm/trunk-jpl/src/m/enum/MatdamageiceEnum.m =================================================================== --- /issm/trunk-jpl/src/m/enum/MatdamageiceEnum.m (revision 17757) +++ /issm/trunk-jpl/src/m/enum/MatdamageiceEnum.m (revision 17757) @@ -0,0 +1,11 @@ +function macro=MatdamageiceEnum() +%MATDAMAGEICEENUM - Enum of Matdamageice +% +% WARNING: DO NOT MODIFY THIS FILE +% this file has been automatically generated by src/c/shared/Enum/Synchronize.sh +% Please read src/c/shared/Enum/README for more information +% +% Usage: +% macro=MatdamageiceEnum() + +macro=StringToEnum('Matdamageice'); Index: /issm/trunk-jpl/src/m/mech/mechanicalproperties.m =================================================================== --- /issm/trunk-jpl/src/m/mech/mechanicalproperties.m (revision 17756) +++ /issm/trunk-jpl/src/m/mech/mechanicalproperties.m (revision 17757) @@ -116,5 +116,5 @@ stress.principalvalue2=valuesstress(:,2); stress.principalaxis2=directionsstress(:,3:4); -stress.effectivevalue=1/sqrt(2)*sqrt(stress.xx.^2+stress.yy.^2+2*stress.xy.^2); % effective shear stress +stress.effectivevalue=1/sqrt(2)*sqrt(stress.xx.^2+stress.yy.^2+2*stress.xy.^2); % effective shear stress (sqrt(J2)) md.results.stress=stress;