Index: /issm/trunk/src/c/EnumDefinitions/EnumDefinitions.cpp =================================================================== --- /issm/trunk/src/c/EnumDefinitions/EnumDefinitions.cpp (revision 2713) +++ /issm/trunk/src/c/EnumDefinitions/EnumDefinitions.cpp (revision 2714) @@ -7,103 +7,104 @@ /*Datasets: */ -int DatasetsEnum(void){ return 100; } -int ElementsEnum(void){ return 101; } -int NodesEnum(void){ return 102; } -int ConstraintsEnum(void){ return 103; } -int LoadsEnum(void){ return 104; } -int MaterialsEnum(void){ return 105; } -int ParametersEnum(void){ return 106; } -int ResultsEnum(void){ return 107; } +int DatasetsEnum(void){ return 100; } +int ElementsEnum(void){ return 101; } +int NodesEnum(void){ return 102; } +int ConstraintsEnum(void){ return 103; } +int LoadsEnum(void){ return 104; } +int MaterialsEnum(void){ return 105; } +int ParametersEnum(void){ return 106; } +int ResultsEnum(void){ return 107; } /*ANALYSIS TYPES: */ -int AnalysisEnum(void){ return 200; } +int AnalysisEnum(void){ return 200; } //diagnostic -int DiagnosticAnalysisEnum(void){ return 210; } -int HorizAnalysisEnum(void){ return 211; } -int StokesAnalysisEnum(void){ return 212; } -int HutterAnalysisEnum(void){ return 213; } -int VertAnalysisEnum(void){ return 214; } +int DiagnosticAnalysisEnum(void){ return 210; } +int HorizAnalysisEnum(void){ return 211; } +int StokesAnalysisEnum(void){ return 212; } +int HutterAnalysisEnum(void){ return 213; } +int VertAnalysisEnum(void){ return 214; } //control -int ControlAnalysisEnum(void){ return 220; } -int AdjointAnalysisEnum(void){ return 221; } -int InverseAnalysisEnum(void){ return 222; } -int GradientAnalysisEnum(void){ return 223; } +int ControlAnalysisEnum(void){ return 220; } +int AdjointAnalysisEnum(void){ return 221; } +int InverseAnalysisEnum(void){ return 222; } +int GradientAnalysisEnum(void){ return 223; } //thermal -int ThermalAnalysisEnum(void){ return 230; } +int ThermalAnalysisEnum(void){ return 230; } //transient -int TransientAnalysisEnum(void){ return 231; } +int TransientAnalysisEnum(void){ return 231; } //slope -int SlopeComputeAnalysisEnum(void){ return 240; } -int SurfaceXAnalysisEnum(void){ return 241; } -int SurfaceYAnalysisEnum(void){ return 242; } -int BedXAnalysisEnum(void){ return 243; } -int BedYAnalysisEnum(void){ return 244; } +int SlopeComputeAnalysisEnum(void){ return 240; } +int SurfaceXAnalysisEnum(void){ return 241; } +int SurfaceYAnalysisEnum(void){ return 242; } +int BedXAnalysisEnum(void){ return 243; } +int BedYAnalysisEnum(void){ return 244; } //prognostic -int PrognosticAnalysisEnum(void){ return 250; } +int PrognosticAnalysisEnum(void){ return 250; } +int BalancedthicknessAnalysisEnum(void){ return 251; } //melting -int MeltingAnalysisEnum(void){ return 260; } +int MeltingAnalysisEnum(void){ return 260; } //mesh2grid -int Mesh2gridAnalysisEnum(void){ return 270; } +int Mesh2gridAnalysisEnum(void){ return 270; } //parameters -int ParametersAnalysisEnum(void){ return 280; } +int ParametersAnalysisEnum(void){ return 280; } //steadystate -int SteadystateAnalysisEnum(void){ return 281; } +int SteadystateAnalysisEnum(void){ return 281; } //none -int NoneAnalysisEnum(void){ return 290; } +int NoneAnalysisEnum(void){ return 290; } /*Formulations: */ -int FormulationEnum(void){ return 300; } -int NoneFormulationEnum(void){ return 301; } -int HutterFormulationEnum(void){ return 302; } -int MacAyealFormulationEnum(void){ return 303; } -int PattynFormulationEnum(void){ return 304; } -int StokesFormulationEnum(void){ return 305; } +int FormulationEnum(void){ return 300; } +int NoneFormulationEnum(void){ return 301; } +int HutterFormulationEnum(void){ return 302; } +int MacAyealFormulationEnum(void){ return 303; } +int PattynFormulationEnum(void){ return 304; } +int StokesFormulationEnum(void){ return 305; } /*OBJECTS: */ -int ObjectEnum(void){ return 400; } +int ObjectEnum(void){ return 400; } /*Elements: */ -int ElementEnum(void){ return 410; } -int TriaEnum(void){ return 411; } -int PentaEnum(void){ return 412; } -int SingEnum(void){ return 414; } -int BeamEnum(void){ return 415; } +int ElementEnum(void){ return 410; } +int TriaEnum(void){ return 411; } +int PentaEnum(void){ return 412; } +int SingEnum(void){ return 414; } +int BeamEnum(void){ return 415; } /*Grids: */ -int NodeEnum(void){ return 420; } +int NodeEnum(void){ return 420; } /*Loads: */ -int LoadEnum(void){ return 430; } -int IcefrontEnum(void){ return 431; } -int RiftfrontEnum(void){ return 432; } -int PenpairEnum(void){ return 433; } -int PengridEnum(void){ return 434; } +int LoadEnum(void){ return 430; } +int IcefrontEnum(void){ return 431; } +int RiftfrontEnum(void){ return 432; } +int PenpairEnum(void){ return 433; } +int PengridEnum(void){ return 434; } /*Materials: */ -int MaterialEnum(void){ return 440; } -int MaticeEnum(void){ return 441; } -int MatparEnum(void){ return 442; } -int NumparEnum(void){ return 443; } +int MaterialEnum(void){ return 440; } +int MaticeEnum(void){ return 441; } +int MatparEnum(void){ return 442; } +int NumparEnum(void){ return 443; } /*Inputs: */ -int InputEnum(void){ return 450; } +int InputEnum(void){ return 450; } /*Params: */ -int ParamEnum(void){ return 460; } +int ParamEnum(void){ return 460; } /*Results: */ -int ResultEnum(void){ return 470; } +int ResultEnum(void){ return 470; } /*Rgb: */ -int RgbEnum(void){ return 480; } +int RgbEnum(void){ return 480; } /*Spc: */ -int SpcEnum(void){ return 490; } +int SpcEnum(void){ return 490; } /*DofVec: */ -int DofVecEnum(void){ return 495; } +int DofVecEnum(void){ return 495; } /*GEOGRAPHY:*/ -int GeographyEnum(void){ return 500; } -int IceSheetEnum(void){ return 502; } -int IceShelfEnum(void){ return 502; } +int GeographyEnum(void){ return 500; } +int IceSheetEnum(void){ return 502; } +int IceShelfEnum(void){ return 502; } /*FILL:*/ -int WaterEnum(void){ return 601; } -int IceEnum(void){ return 602; } -int AirEnum(void){ return 603; } -int MelangeEnum(void){ return 604; } +int WaterEnum(void){ return 601; } +int IceEnum(void){ return 602; } +int AirEnum(void){ return 603; } +int MelangeEnum(void){ return 604; } /*functions on enums: */ Index: /issm/trunk/src/c/EnumDefinitions/EnumDefinitions.h =================================================================== --- /issm/trunk/src/c/EnumDefinitions/EnumDefinitions.h (revision 2713) +++ /issm/trunk/src/c/EnumDefinitions/EnumDefinitions.h (revision 2714) @@ -42,4 +42,5 @@ //prognostic int PrognosticAnalysisEnum(void); +int BalancedthicknessAnalysisEnum(void); //melting int MeltingAnalysisEnum(void); Index: /issm/trunk/src/c/Makefile.am =================================================================== --- /issm/trunk/src/c/Makefile.am (revision 2713) +++ /issm/trunk/src/c/Makefile.am (revision 2714) @@ -218,4 +218,8 @@ ./ModelProcessorx/Prognostic/CreateLoadsPrognostic.cpp\ ./ModelProcessorx/Prognostic/CreateParametersPrognostic.cpp\ + ./ModelProcessorx/Balancedthickness/CreateElementsNodesAndMaterialsBalancedthickness.cpp\ + ./ModelProcessorx/Balancedthickness/CreateConstraintsBalancedthickness.cpp\ + ./ModelProcessorx/Balancedthickness/CreateLoadsBalancedthickness.cpp\ + ./ModelProcessorx/Balancedthickness/CreateParametersBalancedthickness.cpp\ ./ModelProcessorx/Qmu/CreateParametersQmu.cpp\ ./Dofx/Dofx.h\ @@ -517,4 +521,8 @@ ./ModelProcessorx/Prognostic/CreateLoadsPrognostic.cpp\ ./ModelProcessorx/Prognostic/CreateParametersPrognostic.cpp\ + ./ModelProcessorx/Balancedthickness/CreateElementsNodesAndMaterialsBalancedthickness.cpp\ + ./ModelProcessorx/Balancedthickness/CreateConstraintsBalancedthickness.cpp\ + ./ModelProcessorx/Balancedthickness/CreateLoadsBalancedthickness.cpp\ + ./ModelProcessorx/Balancedthickness/CreateParametersBalancedthickness.cpp\ ./ModelProcessorx/Qmu/CreateParametersQmu.cpp\ ./Dofx/Dofx.h\ @@ -616,4 +624,5 @@ ./parallel/ProcessResults.cpp\ ./parallel/prognostic_core.cpp\ + ./parallel/balancedthickness_core.cpp\ ./parallel/transient_core.cpp\ ./parallel/transient_core_2d.cpp\ @@ -627,5 +636,5 @@ bin_PROGRAMS = else -bin_PROGRAMS = diagnostic.exe thermal.exe prognostic.exe transient.exe steadystate.exe +bin_PROGRAMS = diagnostic.exe thermal.exe prognostic.exe balancedthickness.exe transient.exe steadystate.exe endif @@ -644,5 +653,7 @@ prognostic_exe_CXXFLAGS= -fPIC -D_PARALLEL_ +balancedthickness_exe_SOURCES = parallel/balancedthickness.cpp +balancedthickness_exe_CXXFLAGS= -fPIC -D_PARALLEL_ + transient_exe_SOURCES = parallel/transient.cpp transient_exe_CXXFLAGS= -fPIC -D_PARALLEL_ - Index: /issm/trunk/src/c/ModelProcessorx/CreateDataSets.cpp =================================================================== --- /issm/trunk/src/c/ModelProcessorx/CreateDataSets.cpp (revision 2713) +++ /issm/trunk/src/c/ModelProcessorx/CreateDataSets.cpp (revision 2714) @@ -86,6 +86,16 @@ } + else if (iomodel->analysis_type==BalancedthicknessAnalysisEnum()){ + + CreateElementsNodesAndMaterialsBalancedthickness(pelements,pnodes,pmaterials, iomodel,iomodel_handle); + CreateConstraintsBalancedthickness(pconstraints,iomodel,iomodel_handle); + CreateLoadsBalancedthickness(ploads,iomodel,iomodel_handle); + CreateParametersBalancedthickness(pparameters,iomodel,iomodel_handle); + + } else{ + throw ErrorException(__FUNCT__,exprintf("%s%i%s%i%s"," analysis_type: ",iomodel->analysis_type," sub_analysis_type: ",iomodel->sub_analysis_type," not supported yet!")); + } CreateParametersQmu(pparameters,iomodel,iomodel_handle); Index: /issm/trunk/src/c/ModelProcessorx/IoModel.h =================================================================== --- /issm/trunk/src/c/ModelProcessorx/IoModel.h (revision 2713) +++ /issm/trunk/src/c/ModelProcessorx/IoModel.h (revision 2714) @@ -247,4 +247,10 @@ void CreateParametersPrognostic(DataSet** pparameters,IoModel* iomodel,ConstDataHandle iomodel_handle); + /*balancedthickness:*/ + void CreateElementsNodesAndMaterialsBalancedthickness(DataSet** pelements,DataSet** pnodes, DataSet** pmaterials, IoModel* iomodel,ConstDataHandle iomodel_handle); + void CreateConstraintsBalancedthickness(DataSet** pconstraints,IoModel* iomodel,ConstDataHandle iomodel_handle); + void CreateLoadsBalancedthickness(DataSet** ploads, IoModel* iomodel, ConstDataHandle iomodel_handle); + void CreateParametersBalancedthickness(DataSet** pparameters,IoModel* iomodel,ConstDataHandle iomodel_handle); + /*qmu: */ void CreateParametersQmu(DataSet** pparameters,IoModel* iomodel,ConstDataHandle iomodel_handle); Index: /issm/trunk/src/c/io/FetchParams.cpp =================================================================== --- /issm/trunk/src/c/io/FetchParams.cpp (revision 2713) +++ /issm/trunk/src/c/io/FetchParams.cpp (revision 2714) @@ -66,5 +66,5 @@ if (M==0 | N==0){ - throw ErrorException(__FUNCT__,exprintf("%s%s%i%s%i%s%i%s",__FUNCT__," error message: array in parameters structure field ",count," is of size (",M,",",N,")")); + throw ErrorException(__FUNCT__,exprintf("%s%s%i (%s) %s%i%s%i%s",__FUNCT__," error message: array in parameters structure field ",count,name," is of size (",M,",",N,")")); } Index: /issm/trunk/src/c/objects/Penta.cpp =================================================================== --- /issm/trunk/src/c/objects/Penta.cpp (revision 2713) +++ /issm/trunk/src/c/objects/Penta.cpp (revision 2714) @@ -321,4 +321,8 @@ CreateKMatrixPrognostic( Kgg,inputs,analysis_type,sub_analysis_type); } + else if (analysis_type==BalancedthicknessAnalysisEnum()){ + + CreateKMatrixBalancedthickness( Kgg,inputs,analysis_type,sub_analysis_type); + } else if (analysis_type==ThermalAnalysisEnum()){ @@ -332,4 +336,27 @@ throw ErrorException(__FUNCT__,exprintf("%s%i%s\n","analysis: ",analysis_type," not supported yet")); } + +} +/*}}}*/ +/*FUNCTION CreateKMatrixBalancedthickness {{{1*/ +#undef __FUNCT__ +#define __FUNCT__ "Penta::CreateKMatrixBalancedthickness" + +void Penta::CreateKMatrixBalancedthickness(Mat Kgg,void* inputs,int analysis_type,int sub_analysis_type){ + + /*Collapsed formulation: */ + Tria* tria=NULL; + + /*If on water, skip: */ + if(onwater)return; + + /*Is this element on the bed? :*/ + if(!onbed)return; + + /*Spawn Tria element from the base of the Penta: */ + tria=(Tria*)SpawnTria(0,1,2); //grids 0, 1 and 2 make the new tria. + tria->CreateKMatrix(Kgg,inputs, analysis_type,sub_analysis_type); + delete tria; + return; } @@ -1379,4 +1406,8 @@ CreatePVectorPrognostic( pg,inputs,analysis_type,sub_analysis_type); } + else if (analysis_type==BalancedthicknessAnalysisEnum()){ + + CreatePVectorBalancedthickness( pg,inputs,analysis_type,sub_analysis_type); + } else if (analysis_type==ThermalAnalysisEnum()){ @@ -1391,4 +1422,26 @@ } +} +/*}}}*/ +/*FUNCTION CreatePVectorBalancedthickness {{{1*/ +#undef __FUNCT__ +#define __FUNCT__ "Penta::CreatePVectorBalancedthickness" + +void Penta::CreatePVectorBalancedthickness( Vec pg, void* inputs, int analysis_type,int sub_analysis_type){ + + /*Collapsed formulation: */ + Tria* tria=NULL; + + /*If on water, skip: */ + if(onwater)return; + + /*Is this element on the bed? :*/ + if(!onbed)return; + + /*Spawn Tria element from the base of the Penta: */ + tria=(Tria*)SpawnTria(0,1,2); //grids 0, 1 and 2 make the new tria. + tria->CreatePVector(pg,inputs, analysis_type,sub_analysis_type); + delete tria; + return; } /*}}}*/ Index: /issm/trunk/src/c/objects/Penta.h =================================================================== --- /issm/trunk/src/c/objects/Penta.h (revision 2713) +++ /issm/trunk/src/c/objects/Penta.h (revision 2714) @@ -116,4 +116,6 @@ void CreateKMatrixPrognostic(Mat Kgg,void* vinputs,int analysis_type,int sub_analysis_type); void CreatePVectorPrognostic( Vec pg, void* vinputs, int analysis_type,int sub_analysis_type); + void CreateKMatrixBalancedthickness(Mat Kgg,void* vinputs,int analysis_type,int sub_analysis_type); + void CreatePVectorBalancedthickness( Vec pg, void* vinputs, int analysis_type,int sub_analysis_type); void CreateKMatrixDiagnosticStokes( Mat Kgg, void* vinputs, int analysis_type,int sub_analysis_type); Index: /issm/trunk/src/c/objects/Tria.cpp =================================================================== --- /issm/trunk/src/c/objects/Tria.cpp (revision 2713) +++ /issm/trunk/src/c/objects/Tria.cpp (revision 2714) @@ -298,7 +298,215 @@ } + else if (analysis_type==BalancedthicknessAnalysisEnum()){ + + CreateKMatrixBalancedthickness( Kgg,inputs,analysis_type,sub_analysis_type); + + } else{ throw ErrorException(__FUNCT__,exprintf("%s%i%s\n","analysis: ",analysis_type," not supported yet")); } + +} +/*}}}*/ +/*FUNCTION CreateKMatrixBalancedthickness {{{1*/ +#undef __FUNCT__ +#define __FUNCT__ "Tria::CreateKMatrixBalancedthickness" +void Tria::CreateKMatrixBalancedthickness(Mat Kgg,void* vinputs,int analysis_type,int sub_analysis_type){ + + /* local declarations */ + int i,j; + + /* node data: */ + const int numgrids=3; + const int NDOF1=1; + const int numdof=NDOF1*numgrids; + double xyz_list[numgrids][3]; + int doflist[numdof]; + int numberofdofspernode; + + /* gaussian points: */ + int num_gauss,ig; + double* first_gauss_area_coord = NULL; + double* second_gauss_area_coord = NULL; + double* third_gauss_area_coord = NULL; + double* gauss_weights = NULL; + double gauss_weight; + double gauss_l1l2l3[3]; + + /* matrices: */ + double L[numgrids]; + double B[2][numgrids]; + double Bprime[2][numgrids]; + double DL[2][2]={0.0}; + double DLprime[2][2]={0.0}; + double DL_scalar; + double Ke_gg[numdof][numdof]={0.0};//local element stiffness matrix + double Ke_gg_gaussian[numdof][numdof]={0.0}; //stiffness matrix evaluated at the gaussian point. + double Ke_gg_thickness1[numdof][numdof]={0.0}; //stiffness matrix evaluated at the gaussian point. + double Ke_gg_thickness2[numdof][numdof]={0.0}; //stiffness matrix evaluated at the gaussian point. + + double Jdettria; + + /*input parameters for structural analysis (diagnostic): */ + double vxvy_list[numgrids][2]={0.0}; + double vx_list[numgrids]={0.0}; + double vy_list[numgrids]={0.0}; + double dvx[2]; + double dvy[2]; + double vx,vy; + double dvxdx,dvydy; + double v_gauss[2]={0.0}; + double K[2][2]={0.0}; + double KDL[2][2]={0.0}; + int dofs[2]={0,1}; + int found=0; + + ParameterInputs* inputs=NULL; + + /*recover pointers: */ + inputs=(ParameterInputs*)vinputs; + + /*recover extra inputs from users, at current convergence iteration: */ + found=inputs->Recover("velocity_average",&vxvy_list[0][0],2,dofs,numgrids,(void**)nodes); + if(!found)throw ErrorException(__FUNCT__," could not find velocity_average in inputs!"); + + for(i=0;iartdiff){ + //Get the Jacobian determinant + gauss_l1l2l3[0]=1.0/3.0; gauss_l1l2l3[1]=1.0/3.0; gauss_l1l2l3[2]=1.0/3.0; + GetJacobianDeterminant2d(&Jdettria, &xyz_list[0][0],gauss_l1l2l3); + + //Build K matrix (artificial diffusivity matrix) + v_gauss[0]=1.0/3.0*(vxvy_list[0][0]+vxvy_list[1][0]+vxvy_list[2][0]); + v_gauss[1]=1.0/3.0*(vxvy_list[0][1]+vxvy_list[1][1]+vxvy_list[2][1]); + + K[0][0]=pow(Jdettria,(double).5)/2.0*fabs(v_gauss[0]); + K[1][1]=pow(Jdettria,(double).5)/2.0*fabs(v_gauss[1]); + } + + /* Get gaussian points and weights (make this a statically initialized list of points? fstd): */ + GaussTria( &num_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &gauss_weights, 2); + + /* Start looping on the number of gaussian points: */ + for (ig=0; igartdiff){ + + /* Compute artificial diffusivity */ + KDL[0][0]=DL_scalar*K[0][0]; + KDL[1][1]=DL_scalar*K[1][1]; + + TripleMultiply( &Bprime[0][0],2,numdof,1, + &KDL[0][0],2,2,0, + &Bprime[0][0],2,numdof,0, + &Ke_gg_gaussian[0][0],0); + + /* Add artificial diffusivity matrix */ + for( i=0; iRecover("accumulation",&accumulation_list[0],1,dofs,numgrids,(void**)nodes); + if(!found)throw ErrorException(__FUNCT__," could not find accumulation in inputs!"); + found=inputs->Recover("melting",&melting_list[0],1,dofs,numgrids,(void**)nodes); + if(!found)throw ErrorException(__FUNCT__," could not find melting in inputs!"); + + /* Get node coordinates and dof list: */ + GetElementNodeData( &xyz_list[0][0], nodes, numgrids); + GetDofList(&doflist[0],&numberofdofspernode); + + /* Get gaussian points and weights (make this a statically initialized list of points? fstd): */ + GaussTria( &num_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &gauss_weights, 2); + + /* Start looping on the number of gaussian points: */ + for (ig=0; igGetFormulation(DiagnosticAnalysisEnum(),HutterAnalysisEnum()); fem_sl=model->GetFormulation(SlopeComputeAnalysisEnum()); - fem_p=model->GetFormulation(PrognosticAnalysisEnum()); + if(analysis_type==PrognosticAnalysisEnum()){ + fem_p=model->GetFormulation(PrognosticAnalysisEnum()); + } + if(analysis_type==BalancedthicknessAnalysisEnum()){ + fem_p=model->GetFormulation(BalancedthicknessAnalysisEnum()); + } fem_t=model->GetFormulation(ThermalAnalysisEnum()); Index: /issm/trunk/src/c/parallel/parallel.h =================================================================== --- /issm/trunk/src/c/parallel/parallel.h (revision 2713) +++ /issm/trunk/src/c/parallel/parallel.h (revision 2714) @@ -17,4 +17,5 @@ void diagnostic_core(DataSet* results,Model* model, ParameterInputs* inputs); void prognostic_core(DataSet* results,Model* model, ParameterInputs* inputs); +void balancedthickness_core(DataSet* results,Model* model, ParameterInputs* inputs); void control_core(DataSet* results,Model* model, ParameterInputs* inputs); Index: /issm/trunk/src/c/parallel/prognostic.cpp =================================================================== --- /issm/trunk/src/c/parallel/prognostic.cpp (revision 2713) +++ /issm/trunk/src/c/parallel/prognostic.cpp (revision 2714) @@ -28,5 +28,4 @@ Model* model=NULL; - Vec h_g=NULL; Vec u_g=NULL; double* u_g_serial=NULL; Index: /issm/trunk/src/c/shared/Dofs/DistributeNumDofs.cpp =================================================================== --- /issm/trunk/src/c/shared/Dofs/DistributeNumDofs.cpp (revision 2713) +++ /issm/trunk/src/c/shared/Dofs/DistributeNumDofs.cpp (revision 2714) @@ -56,4 +56,7 @@ numdofs=1; } + else if (analysis_type==BalancedthicknessAnalysisEnum()){ + numdofs=1; + } else throw ErrorException(__FUNCT__,exprintf("%s%i%s"," analysis type: ",analysis_type," not implemented yet!")); Index: /issm/trunk/src/m/classes/public/ismodelselfconsistent.m =================================================================== --- /issm/trunk/src/m/classes/public/ismodelselfconsistent.m (revision 2713) +++ /issm/trunk/src/m/classes/public/ismodelselfconsistent.m (revision 2714) @@ -310,8 +310,22 @@ checknan(md,fields); - %INITIAL TEMPERATURE, MELTING AND ACCUMULATION + %INITIAL TEMPERATURE fields={'temperature','spctemperature(:,2)','observed_temperature'}; checkgreater(md,fields,0) +end + +%BALANCEDTHICKNESS +if md.analysis_type==BalancedthicknessAnalysisEnum + + %VELOCITIES MELTING AND ACCUMULATION + fields={'vx','vy','accumulation','melting'}; + checksize(md,fields,[md.numberofgrids 1]); + checknan(md,fields); + + %SPC + if any(md.spcthickness(find(md.gridonboundary))~=1), + error(['model not consistent: model ' md.name ' should have all the nodes on boundary constrained in field spcthickness']); + end end Index: /issm/trunk/src/m/classes/public/process_solve_options.m =================================================================== --- /issm/trunk/src/m/classes/public/process_solve_options.m (revision 2713) +++ /issm/trunk/src/m/classes/public/process_solve_options.m (revision 2714) @@ -17,5 +17,5 @@ %check solution type is supported -if ~ismemberi(analysis_type,{'diagnostic','prognostic','thermal','steadystate','parameters','transient'}), +if ~ismemberi(analysis_type,{'diagnostic','prognostic','thermal','steadystate','parameters','transient','balancedthickness'}), error(['process_solve_options error message: analysis_type ' analysis_type ' not supported yet!']); else Index: /issm/trunk/src/m/classes/public/solve.m =================================================================== --- /issm/trunk/src/m/classes/public/solve.m (revision 2713) +++ /issm/trunk/src/m/classes/public/solve.m (revision 2714) @@ -71,4 +71,6 @@ md=steadystate(md); +elseif md.analysis_type==BalancedthicknessAnalysisEnum, + md=balancedthickness(md); else error('solution type not supported for this model configuration.'); Index: /issm/trunk/src/m/enum/AnalysisTypeFromEnum.m =================================================================== --- /issm/trunk/src/m/enum/AnalysisTypeFromEnum.m (revision 2713) +++ /issm/trunk/src/m/enum/AnalysisTypeFromEnum.m (revision 2714) @@ -85,4 +85,8 @@ end +if enum==BalancedthicknessAnalysisEnum(), + string='balancedthickness'; +end + if enum==MeltingAnalysisEnum(), string='melting'; Index: /issm/trunk/src/pro/bytscl.m =================================================================== --- /issm/trunk/src/pro/bytscl.m (revision 2713) +++ /issm/trunk/src/pro/bytscl.m (revision 2714) @@ -3,7 +3,8 @@ %Equivalent of the bytscl idl routine. -Min=min(min(value)); -Max=max(max(value)); +Min=min(value(:)); +Max=max(value(:)); Top=255; -value=(Top + 0.9999)*(value - Min)/(Max - Min); +pow=1; +value=(Top + 0.9999)*(value - Min).^(pow)/(Max - Min)^(pow);