Index: /issm/trunk-jpl/src/c/analyses/Balancethickness2Analysis.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/Balancethickness2Analysis.cpp (revision 18550) +++ /issm/trunk-jpl/src/c/analyses/Balancethickness2Analysis.cpp (revision 18551) @@ -34,4 +34,5 @@ iomodel->FetchDataToInput(elements,SurfaceforcingsMassBalanceEnum); iomodel->FetchDataToInput(elements,BalancethicknessThickeningRateEnum); + iomodel->FetchDataToInput(elements,BalancethicknessCmuEnum); }/*}}}*/ void Balancethickness2Analysis::CreateNodes(Nodes* nodes,IoModel* iomodel){/*{{{*/ @@ -144,5 +145,5 @@ for(int i=0;ivalues[i]+=Jdet*gauss->weight*( (ms-mb-dhdt)*basis[i] - -D*db[0]*dbasis[0*numnodes+i] - D*db[1]*dbasis[1*numnodes+i] + //-D*db[0]*dbasis[0*numnodes+i] - D*db[1]*dbasis[1*numnodes+i] ); } @@ -155,5 +156,5 @@ }/*}}}*/ void Balancethickness2Analysis::GetSolutionFromInputs(Vector* solution,Element* element){/*{{{*/ - _error_("not implemented yet"); + element->GetSolutionFromInputsOneDof(solution,SurfaceEnum); }/*}}}*/ void Balancethickness2Analysis::GradientJ(Vector* gradient,Element* element,int control_type,int control_index){/*{{{*/ @@ -161,5 +162,55 @@ }/*}}}*/ void Balancethickness2Analysis::InputUpdateFromSolution(IssmDouble* solution,Element* element){/*{{{*/ - element->InputUpdateFromSolutionOneDof(solution,ThicknessEnum); + + /*Intermediaries*/ + IssmDouble ds[2],s,b,D; + IssmDouble* xyz_list = NULL; + + //element->InputUpdateFromSolutionOneDof(solution,ThicknessEnum); + element->InputUpdateFromSolutionOneDof(solution,SurfaceEnum); + + /*Fetch number of vertices and allocate velocity vectors*/ + int numvertices = element->GetNumberOfVertices(); + IssmDouble* vel_list = xNew(numvertices); + IssmDouble* vx_list = xNew(numvertices); + IssmDouble* vy_list = xNew(numvertices); + + /*Retrieve all inputs and parameters*/ + element->GetVerticesCoordinates(&xyz_list); + Input* D_input = element->GetInput(BalancethicknessDiffusionCoefficientEnum); + Input* H_input = element->GetInput(ThicknessEnum); _assert_(H_input); + Input* s_input = element->GetInput(SurfaceEnum); _assert_(s_input); + Input* b_input = element->GetInput(BaseEnum); _assert_(b_input); + + /*Calculate velocities*/ + Gauss* gauss=element->NewGauss(); + for(int iv=0;ivGaussVertex(iv); + + if(D_input){ + D_input->GetInputValue(&D,gauss); + } + else{ + D = 0.; + } + b_input->GetInputValue(&b,gauss); + s_input->GetInputValue(&s,gauss); + s_input->GetInputDerivativeValue(&ds[0],xyz_list,gauss); + + vx_list[iv] = -1./(s-b)*D*ds[0]; + vy_list[iv] = -1./(s-b)*D*ds[1]; + vel_list[iv] = sqrt(pow(vx_list[iv],2) + pow(vy_list[iv],2)); + } + + /*Add vx and vy as inputs to the tria element: */ + element->AddInput(VxEnum,vx_list,P1Enum); + element->AddInput(VyEnum,vy_list,P1Enum); + element->AddInput(VelEnum,vel_list,P1Enum); + + /*Free ressources:*/ + xDelete(vy_list); + xDelete(vx_list); + xDelete(vel_list); + xDelete(xyz_list); }/*}}}*/ void Balancethickness2Analysis::UpdateConstraints(FemModel* femmodel){/*{{{*/ @@ -172,5 +223,5 @@ /*Intermediaries */ - IssmDouble Gamma,h,mu0,ds[2],Cmu,B,k; + IssmDouble Gamma,h,mu0,ds[2],Cmu,B,k,s,b; const int n = 3; const IssmDouble Hstar = 500.; @@ -185,8 +236,9 @@ /*retrieve what we need: */ element->GetVerticesCoordinates(&xyz_list); - Input* thickness_input = element->GetInput(ThicknessEnum); _assert_(thickness_input); + Input* bed_input = element->GetInput(BaseEnum); _assert_(bed_input); Input* surface_input = element->GetInput(SurfaceEnum); _assert_(surface_input); Input* B_input = element->GetInput(MaterialsRheologyBEnum); _assert_(B_input); Input* k_input = element->GetInput(FrictionCoefficientEnum);_assert_(k_input); + Input* Cmu_input = element->GetInput(BalancethicknessCmuEnum);_assert_(Cmu_input); IssmDouble rhog = element->GetMaterialParameter(MaterialsRhoIceEnum)*element->GetMaterialParameter(ConstantsGEnum); @@ -198,12 +250,16 @@ B_input->GetInputValue(&B,gauss); k_input->GetInputValue(&k,gauss); - thickness_input->GetInputValue(&h,gauss); + //thickness_input->GetInputValue(&h,gauss); surface_input->GetInputDerivativeValue(&ds[0],xyz_list,gauss); + surface_input->GetInputValue(&s,gauss); + bed_input->GetInputValue(&b,gauss); + h = s-b; mu0 = pow(2.,(1-3*n)/(2.*n)) * B; Gamma = pow(rhog,n) * pow(Hstar,2*(n+1)) * pow(Hstar/Lstar,2*(n+1)) * 1./pow(mu0,n); - Cmu = mu0/Hstar/(k*k); - - D[i] = Gamma*pow(ds[0]*ds[0]+ds[1]*ds[1],(n-1)/2)*pow(h,n+1)*(1./(n+2.)*h - Cmu); + //Cmu = pow(mu0,n)/Hstar/(k*k); + Cmu_input->GetInputValue(&Cmu,gauss); + + D[i] = Gamma*pow(ds[0]*ds[0]+ds[1]*ds[1],(n-1)/2)*pow(h,n+1)*(1./(n+2.)*h + Cmu); } Index: /issm/trunk-jpl/src/c/cores/balancethickness2_core.cpp =================================================================== --- /issm/trunk-jpl/src/c/cores/balancethickness2_core.cpp (revision 18550) +++ /issm/trunk-jpl/src/c/cores/balancethickness2_core.cpp (revision 18551) @@ -22,10 +22,13 @@ if(VerboseSolution()) _printf0_("call computational core:\n"); - solutionsequence_linear(femmodel); + //solutionsequence_linear(femmodel); + solutionsequence_nonlinear(femmodel,false); if(save_results){ if(VerboseSolution()) _printf0_(" saving results\n"); - int outputs[1] = {ThicknessEnum}; - femmodel->RequestedOutputsx(&femmodel->results,&outputs[0],1); + const int numoutputs = 5; + //int outputs[numoutputs] = {ThicknessEnum}; + int outputs[numoutputs] = {SurfaceEnum,VxEnum,VyEnum,VelEnum,BalancethicknessDiffusionCoefficientEnum}; + femmodel->RequestedOutputsx(&femmodel->results,&outputs[0],numoutputs); } Index: /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h (revision 18550) +++ /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h (revision 18551) @@ -305,4 +305,5 @@ Balancethickness2MisfitEnum, BalancethicknessDiffusionCoefficientEnum, + BalancethicknessCmuEnum, /*}}}*/ /*Surfaceforcings{{{*/ Index: /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp (revision 18550) +++ /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp (revision 18551) @@ -313,4 +313,5 @@ case Balancethickness2MisfitEnum : return "Balancethickness2Misfit"; case BalancethicknessDiffusionCoefficientEnum : return "BalancethicknessDiffusionCoefficient"; + case BalancethicknessCmuEnum : return "BalancethicknessCmu"; case SurfaceforcingsEnum : return "Surfaceforcings"; case SMBEnum : return "SMB"; Index: /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp (revision 18550) +++ /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp (revision 18551) @@ -319,4 +319,5 @@ else if (strcmp(name,"Balancethickness2Misfit")==0) return Balancethickness2MisfitEnum; else if (strcmp(name,"BalancethicknessDiffusionCoefficient")==0) return BalancethicknessDiffusionCoefficientEnum; + else if (strcmp(name,"BalancethicknessCmu")==0) return BalancethicknessCmuEnum; else if (strcmp(name,"Surfaceforcings")==0) return SurfaceforcingsEnum; else if (strcmp(name,"SMB")==0) return SMBEnum; @@ -382,9 +383,9 @@ else if (strcmp(name,"MeltingAnalysis")==0) return MeltingAnalysisEnum; else if (strcmp(name,"MasstransportAnalysis")==0) return MasstransportAnalysisEnum; - else if (strcmp(name,"MasstransportSolution")==0) return MasstransportSolutionEnum; else stage=4; } if(stage==4){ - if (strcmp(name,"FreeSurfaceBaseAnalysis")==0) return FreeSurfaceBaseAnalysisEnum; + if (strcmp(name,"MasstransportSolution")==0) return MasstransportSolutionEnum; + else if (strcmp(name,"FreeSurfaceBaseAnalysis")==0) return FreeSurfaceBaseAnalysisEnum; else if (strcmp(name,"FreeSurfaceTopAnalysis")==0) return FreeSurfaceTopAnalysisEnum; else if (strcmp(name,"SurfaceNormalVelocity")==0) return SurfaceNormalVelocityEnum; @@ -505,9 +506,9 @@ else if (strcmp(name,"FractionIncrement")==0) return FractionIncrementEnum; else if (strcmp(name,"Friction")==0) return FrictionEnum; - else if (strcmp(name,"Internal")==0) return InternalEnum; else stage=5; } if(stage==5){ - if (strcmp(name,"MassFlux")==0) return MassFluxEnum; + if (strcmp(name,"Internal")==0) return InternalEnum; + else if (strcmp(name,"MassFlux")==0) return MassFluxEnum; else if (strcmp(name,"MeltingOffset")==0) return MeltingOffsetEnum; else if (strcmp(name,"Misfit")==0) return MisfitEnum; @@ -628,9 +629,9 @@ else if (strcmp(name,"MisfitObservation")==0) return MisfitObservationEnum; else if (strcmp(name,"MisfitObservationEnum")==0) return MisfitObservationEnumEnum; - else if (strcmp(name,"MisfitTimeinterpolation")==0) return MisfitTimeinterpolationEnum; else stage=6; } if(stage==6){ - if (strcmp(name,"MisfitWeights")==0) return MisfitWeightsEnum; + if (strcmp(name,"MisfitTimeinterpolation")==0) return MisfitTimeinterpolationEnum; + else if (strcmp(name,"MisfitWeights")==0) return MisfitWeightsEnum; else if (strcmp(name,"MisfitWeightsEnum")==0) return MisfitWeightsEnumEnum; else if (strcmp(name,"SurfaceObservation")==0) return SurfaceObservationEnum; @@ -751,9 +752,9 @@ else if (strcmp(name,"MaterialsYoungModulus")==0) return MaterialsYoungModulusEnum; else if (strcmp(name,"MaterialsTimeRelaxationStress")==0) return MaterialsTimeRelaxationStressEnum; - else if (strcmp(name,"MaterialsTimeRelaxationDamage")==0) return MaterialsTimeRelaxationDamageEnum; else stage=7; } if(stage==7){ - if (strcmp(name,"MaterialsRidgingExponent")==0) return MaterialsRidgingExponentEnum; + if (strcmp(name,"MaterialsTimeRelaxationDamage")==0) return MaterialsTimeRelaxationDamageEnum; + else if (strcmp(name,"MaterialsRidgingExponent")==0) return MaterialsRidgingExponentEnum; else if (strcmp(name,"MaterialsCohesion")==0) return MaterialsCohesionEnum; else if (strcmp(name,"MaterialsInternalFrictionCoef")==0) return MaterialsInternalFrictionCoefEnum; Index: /issm/trunk-jpl/src/m/classes/balancethickness.m =================================================================== --- /issm/trunk-jpl/src/m/classes/balancethickness.m (revision 18550) +++ /issm/trunk-jpl/src/m/classes/balancethickness.m (revision 18551) @@ -9,20 +9,8 @@ thickening_rate = NaN; stabilization = 0; + + Cmu = NaN; end methods - function createxml(obj,fid) % {{{ - fprintf(fid, '\n\n'); - fprintf(fid, '%s\n', ''); - - fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n', '', '
',' thickness constraints (NaN means no constraint) [m] ',''); - fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n', '', '
',' ice thickening rate used in the mass conservation (dh/dt) [m/yr] ',''); - - % balancethickness drop-down (1,2, or 3) - fprintf(fid,'%s\n%s\n%s\n', '','
',' 0: None, 1: SU, 2: SSAs artificial diffusivity, 3:DG '); - fprintf(fid,'%s\n',' '); - fprintf(fid,'%s\n',' '); - fprintf(fid,'%s\n%s\n',' ',''); - - end % }}} function obj = balancethickness(varargin) % {{{ switch nargin @@ -46,4 +34,6 @@ md = checkfield(md,'fieldname','balancethickness.thickening_rate','size',[md.mesh.numberofvertices 1],'NaN',1); md = checkfield(md,'fieldname','balancethickness.stabilization','size',[1 1],'values',[0 1 2 3]); + + md = checkfield(md,'fieldname','balancethickness.Cmu','size',[md.mesh.numberofvertices 1],'NaN',1,'>=',0); end % }}} function disp(obj) % {{{ @@ -62,4 +52,6 @@ WriteData(fid,'object',obj,'fieldname','thickening_rate','format','DoubleMat','mattype',1,'scale',1./yts); WriteData(fid,'object',obj,'fieldname','stabilization','format','Integer'); + + WriteData(fid,'object',obj,'fieldname','Cmu','format','DoubleMat','mattype',1); end % }}} end Index: /issm/trunk-jpl/src/m/enum/BalancethicknessCmuEnum.m =================================================================== --- /issm/trunk-jpl/src/m/enum/BalancethicknessCmuEnum.m (revision 18551) +++ /issm/trunk-jpl/src/m/enum/BalancethicknessCmuEnum.m (revision 18551) @@ -0,0 +1,11 @@ +function macro=BalancethicknessCmuEnum() +%BALANCETHICKNESSCMUENUM - Enum of BalancethicknessCmu +% +% 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=BalancethicknessCmuEnum() + +macro=StringToEnum('BalancethicknessCmu'); Index: /issm/trunk-jpl/src/m/enum/EnumDefinitions.py =================================================================== --- /issm/trunk-jpl/src/m/enum/EnumDefinitions.py (revision 18550) +++ /issm/trunk-jpl/src/m/enum/EnumDefinitions.py (revision 18551) @@ -305,4 +305,5 @@ def Balancethickness2MisfitEnum(): return StringToEnum("Balancethickness2Misfit")[0] def BalancethicknessDiffusionCoefficientEnum(): return StringToEnum("BalancethicknessDiffusionCoefficient")[0] +def BalancethicknessCmuEnum(): return StringToEnum("BalancethicknessCmu")[0] def SurfaceforcingsEnum(): return StringToEnum("Surfaceforcings")[0] def SMBEnum(): return StringToEnum("SMB")[0] Index: /issm/trunk-jpl/src/m/plot/googlemaps.m =================================================================== --- /issm/trunk-jpl/src/m/plot/googlemaps.m (revision 18550) +++ /issm/trunk-jpl/src/m/plot/googlemaps.m (revision 18551) @@ -54,4 +54,6 @@ elseif strcmpi(md.mesh.hemisphere,'s'), EPSGlocal = 'EPSG:3031'; % UPS Antarctica http://www.spatialreference.org/ref/epsg/3031/ +elseif strcmpi(md.mesh.hemisphere,'UTM18N'), + EPSGlocal = 'EPSG:32218'; % UTM 18 N Patagonia else error('field hemisphere should either be ''n'' or ''s'''); Index: /issm/trunk-jpl/src/m/plot/plot_googlemaps.m =================================================================== --- /issm/trunk-jpl/src/m/plot/plot_googlemaps.m (revision 18550) +++ /issm/trunk-jpl/src/m/plot/plot_googlemaps.m (revision 18551) @@ -36,4 +36,6 @@ -1,0,71); else + latlist = md.mesh.lat; %That might work? + lonlist = md.mesh.long; error('field hemisphere should either be ''n'' or ''s'''); end