Index: /issm/trunk-jpl/src/c/classes/AdaptiveMeshRefinement.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/AdaptiveMeshRefinement.cpp (revision 22293) +++ /issm/trunk-jpl/src/c/classes/AdaptiveMeshRefinement.cpp (revision 22294) @@ -112,4 +112,7 @@ if(this->groundingline_distance>0 && !gl_distance) _error_("gl_distance is NULL!\n"); if(this->icefront_distance>0 && !if_distance) _error_("if_distance is NULL!\n"); + /*Attributes verifications*/ + if(this->deviatoricerror_threshold>0 && this->deviatoricerror_groupthresholdthicknesserror_threshold>0 && this->thicknesserror_groupthresholdlag*gl_distance_hmax*std::pow(this->gradation,this->level_max-gmesh->Element(i)->Level()); if_distance_h = this->lag*if_distance_hmax*std::pow(this->gradation,this->level_max-gmesh->Element(i)->Level()); - d_maxerror = 0.01*this->deviatoricerror_maximum;//itapopo definir melhor 0.05 - t_maxerror = 0.01*this->thicknesserror_maximum;//itapopo definir melhor 0.05 + d_maxerror = this->deviatoricerror_groupthreshold*this->deviatoricerror_maximum; + t_maxerror = this->thicknesserror_groupthreshold*this->thicknesserror_maximum; /*Get the sons of the father (sibilings)*/ sons.clear(); Index: /issm/trunk-jpl/src/c/classes/AdaptiveMeshRefinement.h =================================================================== --- /issm/trunk-jpl/src/c/classes/AdaptiveMeshRefinement.h (revision 22293) +++ /issm/trunk-jpl/src/c/classes/AdaptiveMeshRefinement.h (revision 22294) @@ -64,6 +64,8 @@ double icefront_distance; //all elements with distance from ice front <= icefront_distance will be refined double thicknesserror_threshold; //if ==0, it will not be used + double thicknesserror_groupthreshold;//group threshold double thicknesserror_maximum; //max value of the error estimator; in general, it is defined in the first time step. Attention with restart double deviatoricerror_threshold; //if ==0, it will not be used + double deviatoricerror_groupthreshold;//group threshold double deviatoricerror_maximum; //max value of the error estimator; in general, it is defined in the first time step. Attention with restart /*}}}*/ Index: /issm/trunk-jpl/src/c/classes/AmrBamg.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/AmrBamg.cpp (revision 22293) +++ /issm/trunk-jpl/src/c/classes/AmrBamg.cpp (revision 22294) @@ -20,47 +20,49 @@ /*These attributes MUST be setup by FemModel*/ - this->fieldenum = -1;//fieldenum_in; - this->keepmetric = -1;//keepmetric_in; - this->groundingline_resolution = -1;//groundingline_resolution_in; - this->groundingline_distance = -1;//groundingline_distance_in; - this->icefront_resolution = -1;//icefront_resolution_in; - this->icefront_distance = -1;//icefront_distance_in; - this->thicknesserror_resolution = -1;//thicknesserror_resolution_in; - this->thicknesserror_threshold = -1;//thicknesserror_threshold_in; - this->thicknesserror_maximum = -1; - this->deviatoricerror_resolution = -1;//deviatoricerror_resolution_in; - this->deviatoricerror_threshold = -1;//deviatoricerror_threshold_in; - this->deviatoricerror_maximum = -1; + this->fieldenum = -1; + this->keepmetric = -1; + this->groundingline_resolution = -1; + this->groundingline_distance = -1; + this->icefront_resolution = -1; + this->icefront_distance = -1; + this->thicknesserror_resolution = -1; + this->thicknesserror_threshold = -1; + this->thicknesserror_groupthreshold = -1; + this->thicknesserror_maximum = -1; + this->deviatoricerror_resolution = -1; + this->deviatoricerror_threshold = -1; + this->deviatoricerror_groupthreshold = -1; + this->deviatoricerror_maximum = -1; /*Geometry and mesh as NULL*/ - this->geometry = NULL; - this->fathermesh = NULL; - this->previousmesh = NULL; + this->geometry = NULL; + this->fathermesh = NULL; + this->previousmesh = NULL; /*Only initialize options for now (same as bamg.m)*/ - this->options = new BamgOpts(); - this->options->anisomax = 10.e30; - this->options->cutoff = 10.e-5; - this->options->coeff = 1; - this->options->errg = 0.1; - this->options->gradation = -1; //MUST be setup by the FemModel - this->options->Hessiantype = 0; - this->options->maxnbv = 1e6; - this->options->maxsubdiv = 10; - this->options->Metrictype = 0; - this->options->nbjacobi = 1; - this->options->nbsmooth = 3; - this->options->omega = 1.8; - this->options->power = 1; - this->options->verbose = 0; - this->options->Crack = 0; - this->options->KeepVertices = 1; /*!!!!! VERY IMPORTANT !!!!! This avoid numerical errors when remeshing*/ - this->options->splitcorners = 1; - this->options->hmin = -1;/*MUST be setup by the FemModel*/ - this->options->hmax = -1;/*MUST be setup by the FemModel*/ - this->options->err = xNew(1); - this->options->err[0] = -1;/*MUST be setup by the FemModel*/ - this->options->errSize[0] = 1; - this->options->errSize[1] = 1; + this->options = new BamgOpts(); + this->options->anisomax = 10.e30; + this->options->cutoff = 10.e-5; + this->options->coeff = 1; + this->options->errg = 0.1; + this->options->gradation = -1; //MUST be setup by the FemModel + this->options->Hessiantype = 0; + this->options->maxnbv = 1e6; + this->options->maxsubdiv = 10; + this->options->Metrictype = 0; + this->options->nbjacobi = 1; + this->options->nbsmooth = 3; + this->options->omega = 1.8; + this->options->power = 1; + this->options->verbose = 0; + this->options->Crack = 0; + this->options->KeepVertices = 1; /*!!!!! VERY IMPORTANT !!!!! This avoid numerical errors when remeshing*/ + this->options->splitcorners = 1; + this->options->hmin = -1;/*MUST be setup by the FemModel*/ + this->options->hmax = -1;/*MUST be setup by the FemModel*/ + this->options->err = xNew(1); + this->options->err[0] = -1;/*MUST be setup by the FemModel*/ + this->options->errSize[0] = 1; + this->options->errSize[1] = 1; } /*}}}*/ Index: /issm/trunk-jpl/src/c/classes/AmrBamg.h =================================================================== --- /issm/trunk-jpl/src/c/classes/AmrBamg.h (revision 22293) +++ /issm/trunk-jpl/src/c/classes/AmrBamg.h (revision 22294) @@ -21,7 +21,9 @@ IssmDouble thicknesserror_resolution; IssmDouble thicknesserror_threshold; + IssmDouble thicknesserror_groupthreshold; IssmDouble thicknesserror_maximum; IssmDouble deviatoricerror_resolution; IssmDouble deviatoricerror_threshold; + IssmDouble deviatoricerror_groupthreshold; IssmDouble deviatoricerror_maximum; Index: /issm/trunk-jpl/src/c/classes/FemModel.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/FemModel.cpp (revision 22293) +++ /issm/trunk-jpl/src/c/classes/FemModel.cpp (revision 22294) @@ -2424,9 +2424,9 @@ int* my_vertices = NULL; int elementswidth = this->GetElementsWidth();//just tria elements in this version - int amrtype; + int amrtype,basalforcing_model; bool isgroundingline; /*Branch to specific amr depending on requested method*/ - parameters->FindParam(&amrtype,AmrTypeEnum); + this->parameters->FindParam(&amrtype,AmrTypeEnum); switch(amrtype){ #if defined(_HAVE_NEOPZ_) && !defined(_HAVE_ADOLC_) @@ -2539,5 +2539,6 @@ /*Insert MISMIP+ bed topography FIXME it could be stay in another place*/ - if(false) this->BedrockFromMismipPlus(); + this->parameters->FindParam(&basalforcing_model,BasalforcingsEnum); + if(basalforcing_model==MismipFloatingMeltRateEnum) this->BedrockFromMismipPlus(); /*Adjust base, thickness and mask grounded ice leve set*/ @@ -2563,4 +2564,6 @@ /*This was used to Misomip project/simulations*/ + if(VerboseSolution())_printf0_(" call Mismip bedrock adjust module\n"); + IssmDouble x,y,bx,by; int numvertices = this->GetElementsWidth(); @@ -2570,24 +2573,17 @@ for(int el=0;elelements->Size();el++){ Element* element=xDynamicCast(this->elements->GetObjectByOffset(el)); - element->GetVerticesCoordinates(&xyz_list); for(int i=0;i-720.) r[i] = bx+by; + x = xyz_list[3*i+0]; + y = xyz_list[3*i+1]; + bx = -150.-728.8*pow(x/300000.,2)+343.91*pow(x/300000.,4)-50.57*pow(x/300000.,6); + by = 500./(1.+exp((-2./4000.)*(y-80000./2.-24000.)))+500./(1.+exp((2./4000.)*(y-80000./2.+24000.))); + r[i] = max(bx+by,-720.); } - /*insert new bedrock*/ element->AddInput(BedEnum,&r[0],P1Enum); - /*Cleanup*/ xDelete(xyz_list); } - /*Delete*/ xDelete(r); @@ -2596,4 +2592,6 @@ void FemModel::AdjustBaseThicknessAndMask(void){/*{{{*/ + if(VerboseSolution())_printf0_(" call adjust base and thickness module\n"); + int numvertices = this->GetElementsWidth(); IssmDouble rho_water,rho_ice,density,base_float; @@ -4743,7 +4741,9 @@ this->parameters->FindParam(&this->amrbamg->thicknesserror_resolution,AmrThicknessErrorResolutionEnum); this->parameters->FindParam(&this->amrbamg->thicknesserror_threshold,AmrThicknessErrorThresholdEnum); + this->parameters->FindParam(&this->amrbamg->thicknesserror_groupthreshold,AmrThicknessErrorGroupThresholdEnum); this->parameters->FindParam(&this->amrbamg->thicknesserror_maximum,AmrThicknessErrorMaximumEnum); this->parameters->FindParam(&this->amrbamg->deviatoricerror_resolution,AmrDeviatoricErrorResolutionEnum); this->parameters->FindParam(&this->amrbamg->deviatoricerror_threshold,AmrDeviatoricErrorThresholdEnum); + this->parameters->FindParam(&this->amrbamg->deviatoricerror_groupthreshold,AmrDeviatoricErrorGroupThresholdEnum); this->parameters->FindParam(&this->amrbamg->deviatoricerror_maximum,AmrDeviatoricErrorMaximumEnum); /*Set BamgOpts*/ @@ -4814,5 +4814,5 @@ IssmDouble* z = NULL; int* index = NULL; - IssmDouble maxerror,threshold,resolution,length; + IssmDouble maxerror,threshold,groupthreshold,resolution,length; IssmDouble L1,L2,L3; int vid,v1,v2,v3; @@ -4822,13 +4822,15 @@ switch(errorestimator_type){ case ThicknessErrorEstimatorEnum: - threshold = this->amrbamg->thicknesserror_threshold; - resolution = this->amrbamg->thicknesserror_resolution; - maxerror = this->amrbamg->thicknesserror_maximum; + threshold = this->amrbamg->thicknesserror_threshold; + groupthreshold = this->amrbamg->thicknesserror_groupthreshold; + resolution = this->amrbamg->thicknesserror_resolution; + maxerror = this->amrbamg->thicknesserror_maximum; this->ThicknessZZErrorEstimator(&error_elements);//error is serial, but the calculation is parallel break; case DeviatoricStressErrorEstimatorEnum: - threshold = this->amrbamg->deviatoricerror_threshold; - resolution = this->amrbamg->deviatoricerror_resolution; - maxerror = this->amrbamg->deviatoricerror_maximum; + threshold = this->amrbamg->deviatoricerror_threshold; + groupthreshold = this->amrbamg->deviatoricerror_groupthreshold; + resolution = this->amrbamg->deviatoricerror_resolution; + maxerror = this->amrbamg->deviatoricerror_maximum; this->ZZErrorEstimator(&error_elements);//error is serial, but the calculation is parallel break; @@ -4836,4 +4838,5 @@ } if(!error_elements) _error_("error_elements is NULL!\n"); + if(groupthreshold0.01*maxerror) refine=true;//itapopo must be better defined + if(error_vertices[vid]>groupthreshold*maxerror) refine=true; } } @@ -5034,6 +5037,8 @@ this->parameters->FindParam(&this->amr->icefront_distance,AmrIceFrontDistanceEnum); this->parameters->FindParam(&this->amr->thicknesserror_threshold,AmrThicknessErrorThresholdEnum); + this->parameters->FindParam(&this->amr->thicknesserror_groupthreshold,AmrThicknessErrorGroupThresholdEnum); this->parameters->FindParam(&this->amr->thicknesserror_maximum,AmrThicknessErrorMaximumEnum); this->parameters->FindParam(&this->amr->deviatoricerror_threshold,AmrDeviatoricErrorThresholdEnum); + this->parameters->FindParam(&this->amr->deviatoricerror_groupthreshold,AmrDeviatoricErrorGroupThresholdEnum); this->parameters->FindParam(&this->amr->deviatoricerror_maximum,AmrDeviatoricErrorMaximumEnum); if(my_rank==0){ Index: /issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateParameters.cpp =================================================================== --- /issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateParameters.cpp (revision 22293) +++ /issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateParameters.cpp (revision 22294) @@ -124,5 +124,16 @@ iomodel->FindConstant(&amr_frequency,"md.transient.amr_frequency"); if(solution_type==TransientSolutionEnum && amr_frequency){ + /*Load common amr parameters*/ parameters->AddObject(iomodel->CopyConstantObject("md.amr.type",AmrTypeEnum)); + parameters->AddObject(iomodel->CopyConstantObject("md.amr.groundingline_distance",AmrGroundingLineDistanceEnum)); + parameters->AddObject(iomodel->CopyConstantObject("md.amr.icefront_distance",AmrIceFrontDistanceEnum)); + parameters->AddObject(iomodel->CopyConstantObject("md.amr.thicknesserror_threshold",AmrThicknessErrorThresholdEnum)); + parameters->AddObject(iomodel->CopyConstantObject("md.amr.thicknesserror_groupthreshold",AmrThicknessErrorGroupThresholdEnum)); + parameters->AddObject(iomodel->CopyConstantObject("md.amr.thicknesserror_maximum",AmrThicknessErrorMaximumEnum)); + parameters->AddObject(iomodel->CopyConstantObject("md.amr.deviatoricerror_threshold",AmrDeviatoricErrorThresholdEnum)); + parameters->AddObject(iomodel->CopyConstantObject("md.amr.deviatoricerror_groupthreshold",AmrDeviatoricErrorGroupThresholdEnum)); + parameters->AddObject(iomodel->CopyConstantObject("md.amr.deviatoricerror_maximum",AmrDeviatoricErrorMaximumEnum)); + parameters->AddObject(iomodel->CopyConstantObject("md.amr.restart",AmrRestartEnum)); + /*Load specific amr parameters*/ iomodel->FindConstant(&amrtype,"md.amr.type"); switch(amrtype){ @@ -132,11 +143,4 @@ parameters->AddObject(iomodel->CopyConstantObject("md.amr.gradation",AmrGradationEnum)); parameters->AddObject(iomodel->CopyConstantObject("md.amr.lag",AmrLagEnum)); - parameters->AddObject(iomodel->CopyConstantObject("md.amr.groundingline_distance",AmrGroundingLineDistanceEnum)); - parameters->AddObject(iomodel->CopyConstantObject("md.amr.icefront_distance",AmrIceFrontDistanceEnum)); - parameters->AddObject(iomodel->CopyConstantObject("md.amr.thicknesserror_threshold",AmrThicknessErrorThresholdEnum)); - parameters->AddObject(iomodel->CopyConstantObject("md.amr.thicknesserror_maximum",AmrThicknessErrorMaximumEnum)); - parameters->AddObject(iomodel->CopyConstantObject("md.amr.deviatoricerror_threshold",AmrDeviatoricErrorThresholdEnum)); - parameters->AddObject(iomodel->CopyConstantObject("md.amr.deviatoricerror_maximum",AmrDeviatoricErrorMaximumEnum)); - parameters->AddObject(iomodel->CopyConstantObject("md.amr.restart",AmrRestartEnum)); break; #endif @@ -150,14 +154,7 @@ parameters->AddObject(iomodel->CopyConstantObject("md.amr.gradation",AmrGradationEnum)); parameters->AddObject(iomodel->CopyConstantObject("md.amr.groundingline_resolution",AmrGroundingLineResolutionEnum)); - parameters->AddObject(iomodel->CopyConstantObject("md.amr.groundingline_distance",AmrGroundingLineDistanceEnum)); parameters->AddObject(iomodel->CopyConstantObject("md.amr.icefront_resolution",AmrIceFrontResolutionEnum)); - parameters->AddObject(iomodel->CopyConstantObject("md.amr.icefront_distance",AmrIceFrontDistanceEnum)); parameters->AddObject(iomodel->CopyConstantObject("md.amr.thicknesserror_resolution",AmrThicknessErrorResolutionEnum)); - parameters->AddObject(iomodel->CopyConstantObject("md.amr.thicknesserror_threshold",AmrThicknessErrorThresholdEnum)); - parameters->AddObject(iomodel->CopyConstantObject("md.amr.thicknesserror_maximum",AmrThicknessErrorMaximumEnum)); parameters->AddObject(iomodel->CopyConstantObject("md.amr.deviatoricerror_resolution",AmrDeviatoricErrorResolutionEnum)); - parameters->AddObject(iomodel->CopyConstantObject("md.amr.deviatoricerror_threshold",AmrDeviatoricErrorThresholdEnum)); - parameters->AddObject(iomodel->CopyConstantObject("md.amr.deviatoricerror_maximum",AmrDeviatoricErrorMaximumEnum)); - parameters->AddObject(iomodel->CopyConstantObject("md.amr.restart",AmrRestartEnum)); /*Convert fieldname to enum and put it in params*/ iomodel->FindConstant(&fieldname,"md.amr.fieldname"); Index: /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h (revision 22293) +++ /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h (revision 22294) @@ -901,7 +901,9 @@ AmrThicknessErrorResolutionEnum, AmrThicknessErrorThresholdEnum, + AmrThicknessErrorGroupThresholdEnum, AmrThicknessErrorMaximumEnum, AmrDeviatoricErrorResolutionEnum, AmrDeviatoricErrorThresholdEnum, + AmrDeviatoricErrorGroupThresholdEnum, AmrDeviatoricErrorMaximumEnum, DeviatoricStressErrorEstimatorEnum, Index: /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp (revision 22293) +++ /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp (revision 22294) @@ -874,7 +874,9 @@ case AmrThicknessErrorResolutionEnum : return "AmrThicknessErrorResolution"; case AmrThicknessErrorThresholdEnum : return "AmrThicknessErrorThreshold"; + case AmrThicknessErrorGroupThresholdEnum : return "AmrThicknessErrorGroupThreshold"; case AmrThicknessErrorMaximumEnum : return "AmrThicknessErrorMaximum"; case AmrDeviatoricErrorResolutionEnum : return "AmrDeviatoricErrorResolution"; case AmrDeviatoricErrorThresholdEnum : return "AmrDeviatoricErrorThreshold"; + case AmrDeviatoricErrorGroupThresholdEnum : return "AmrDeviatoricErrorGroupThreshold"; case AmrDeviatoricErrorMaximumEnum : return "AmrDeviatoricErrorMaximum"; case DeviatoricStressErrorEstimatorEnum : return "DeviatoricStressErrorEstimator"; Index: /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp (revision 22293) +++ /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp (revision 22294) @@ -895,7 +895,9 @@ else if (strcmp(name,"AmrThicknessErrorResolution")==0) return AmrThicknessErrorResolutionEnum; else if (strcmp(name,"AmrThicknessErrorThreshold")==0) return AmrThicknessErrorThresholdEnum; + else if (strcmp(name,"AmrThicknessErrorGroupThreshold")==0) return AmrThicknessErrorGroupThresholdEnum; else if (strcmp(name,"AmrThicknessErrorMaximum")==0) return AmrThicknessErrorMaximumEnum; else if (strcmp(name,"AmrDeviatoricErrorResolution")==0) return AmrDeviatoricErrorResolutionEnum; else if (strcmp(name,"AmrDeviatoricErrorThreshold")==0) return AmrDeviatoricErrorThresholdEnum; + else if (strcmp(name,"AmrDeviatoricErrorGroupThreshold")==0) return AmrDeviatoricErrorGroupThresholdEnum; else if (strcmp(name,"AmrDeviatoricErrorMaximum")==0) return AmrDeviatoricErrorMaximumEnum; else if (strcmp(name,"DeviatoricStressErrorEstimator")==0) return DeviatoricStressErrorEstimatorEnum; @@ -996,10 +998,10 @@ else if (strcmp(name,"HydrologySolution")==0) return HydrologySolutionEnum; 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=9; } if(stage==9){ - if (strcmp(name,"FreeSurfaceBaseAnalysis")==0) return FreeSurfaceBaseAnalysisEnum; + if (strcmp(name,"MasstransportAnalysis")==0) return MasstransportAnalysisEnum; + else 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,"ExtrudeFromBaseAnalysis")==0) return ExtrudeFromBaseAnalysisEnum; Index: /issm/trunk-jpl/src/m/classes/amr.js =================================================================== --- /issm/trunk-jpl/src/m/classes/amr.js (revision 22293) +++ /issm/trunk-jpl/src/m/classes/amr.js (revision 22294) @@ -7,20 +7,22 @@ //methods this.setdefaultparameters = function(){// {{{ - this.hmin = 100.; - this.hmax = 100.e3; - this.fieldname = "Vel"; - this.err = 3.; - this.keepmetric = 1; - this.gradation = 1.5; - this.groundingline_resolution = 500.; - this.groundingline_distance = 0; - this.icefront_resolution = 500; - this.icefront_distance = 0; - this.thicknesserror_resolution = 500; - this.thicknesserror_threshold = 0; - this.thicknesserror_maximum = 0; - this.deviatoricerror_resolution = 500; - this.deviatoricerror_threshold = 0; - this.deviatoricerror_maximum = 0; + this.hmin = 100.; + this.hmax = 100.e3; + this.fieldname = "Vel"; + this.err = 3.; + this.keepmetric = 1; + this.gradation = 1.5; + this.groundingline_resolution = 500.; + this.groundingline_distance = 0; + this.icefront_resolution = 500; + this.icefront_distance = 0; + this.thicknesserror_resolution = 500; + this.thicknesserror_threshold = 0; + this.thicknesserror_groupthreshold = 0; + this.thicknesserror_maximum = 0; + this.deviatoricerror_resolution = 500; + this.deviatoricerror_threshold = 0; + this.deviatoricerror_groupthreshold = 0; + this.deviatoricerror_maximum = 0; }// }}} this.disp= function(){// {{{ @@ -37,7 +39,9 @@ fielddisplay(this,'thicknesserror_resolution','element length when thickness error estimator is used'); fielddisplay(this,'thicknesserror_threshold','maximum threshold thickness error permitted'); + fielddisplay(this,'thicknesserror_groupthreshold','maximum group threshold thickness error permitted'); fielddisplay(this,'thicknesserror_maximum','maximum thickness error permitted'); fielddisplay(this,'deviatoricerror_resolution','element length when deviatoric stress error estimator is used'); fielddisplay(this,'deviatoricerror_threshold','maximum threshold deviatoricstress error permitted'); + fielddisplay(this,'deviatoricerror_groupthreshold','maximum group threshold deviatoricstress error permitted'); fielddisplay(this,'deviatoricerror_maximum','maximum deviatoricstress error permitted'); }// }}} @@ -57,7 +61,9 @@ checkfield(md,'fieldname','amr.thicknesserror_resolution','numel',[1],'>',0,'<',this.hmax,'NaN',1); checkfield(md,'fieldname','amr.thicknesserror_threshold','numel',[1],'>=',0,'<=',1,'NaN',1); + checkfield(md,'fieldname','amr.thicknesserror_groupthreshold','numel',[1],'>=',0,'<=',1,'NaN',1); checkfield(md,'fieldname','amr.thicknesserror_maximum','numel',[1],'>=',0,'NaN',1,'Inf',1); checkfield(md,'fieldname','amr.deviatoricerror_resolution','numel',[1],'>',0,'<',this.hmax,'NaN',1); checkfield(md,'fieldname','amr.deviatoricerror_threshold','numel',[1],'>=',0,'<=',1,'NaN',1); + checkfield(md,'fieldname','amr.deviatoricerror_groupthreshold','numel',[1],'>=',0,'<=',1,'NaN',1); checkfield(md,'fieldname','amr.deviatoricerror_maximum','numel',[1],'>=',0,'NaN',1,'Inf',1); } // }}} @@ -76,7 +82,9 @@ WriteData(fid,prefix,'object',this,'fieldname','thicknesserror_resolution','format','Double'); WriteData(fid,prefix,'object',this,'fieldname','thicknesserror_threshold','format','Double'); + WriteData(fid,prefix,'object',this,'fieldname','thicknesserror_groupthreshold','format','Double'); WriteData(fid,prefix,'object',this,'fieldname','thicknesserror_maximum','format','Double'); WriteData(fid,prefix,'object',this,'fieldname','deviatoricerror_resolution','format','Double'); - WriteData(fid,prefix,'object',this,'fieldname','deviatoricerror_threshold','format','Double'); + writedata(fid,prefix,'object',this,'fieldname','deviatoricerror_threshold','format','double'); + writedata(fid,prefix,'object',this,'fieldname','deviatoricerror_groupthreshold','format','double'); WriteData(fid,prefix,'object',this,'fieldname','deviatoricerror_maximum','format','Double'); }//}}} @@ -85,20 +93,22 @@ //properties // {{{ - this.hmin = 0.; - this.hmax = 0.; - this.fieldname = ""; - this.err = 0.; - this.keepmetric = 0; - this.gradation = 0.; - this.groundingline_resolution = 0.; - this.groundingline_distance = 0.; - this.icefront_resolution = 0.; - this.icefront_distance = 0.; - this.thicknesserror_resolution = 0.; - this.thicknesserror_threshold = 0.; - this.thicknesserror_maximum = 0.; - this.deviatoricerror_resolution = 0.; - this.deviatoricerror_threshold = 0.; - this.deviatoricerror_maximum = 0.; + this.hmin = 0.; + this.hmax = 0.; + this.fieldname = ""; + this.err = 0.; + this.keepmetric = 0; + this.gradation = 0.; + this.groundingline_resolution = 0.; + this.groundingline_distance = 0.; + this.icefront_resolution = 0.; + this.icefront_distance = 0.; + this.thicknesserror_resolution = 0.; + this.thicknesserror_threshold = 0.; + this.thicknesserror_groupthreshold = 0.; + this.thicknesserror_maximum = 0.; + this.deviatoricerror_resolution = 0.; + this.deviatoricerror_threshold = 0.; + this.deviatoricerror_groupthreshold = 0.; + this.deviatoricerror_maximum = 0.; this.setdefaultparameters(); Index: /issm/trunk-jpl/src/m/classes/amr.m =================================================================== --- /issm/trunk-jpl/src/m/classes/amr.m (revision 22293) +++ /issm/trunk-jpl/src/m/classes/amr.m (revision 22294) @@ -18,7 +18,9 @@ thicknesserror_resolution = 0.; thicknesserror_threshold = 0.; + thicknesserror_groupthreshold = 0.; thicknesserror_maximum = 0.; deviatoricerror_resolution = 0.; deviatoricerror_threshold = 0.; + deviatoricerror_groupthreshold = 0.; deviatoricerror_maximum = 0.; restart=0.; @@ -89,7 +91,9 @@ self.thicknesserror_resolution=500.; self.thicknesserror_threshold=0.; + self.thicknesserror_groupthreshold=0.; self.thicknesserror_maximum=0.; self.deviatoricerror_resolution=500.; self.deviatoricerror_threshold=0.; + self.deviatoricerror_groupthreshold=0.; self.deviatoricerror_maximum=0.; @@ -111,7 +115,9 @@ md = checkfield(md,'fieldname','amr.thicknesserror_resolution','numel',[1],'>',0,'<',self.hmax,'NaN',1); md = checkfield(md,'fieldname','amr.thicknesserror_threshold','numel',[1],'>=',0,'<=',1,'NaN',1); + md = checkfield(md,'fieldname','amr.thicknesserror_groupthreshold','numel',[1],'>=',0,'<=',1,'NaN',1); md = checkfield(md,'fieldname','amr.thicknesserror_maximum','numel',[1],'>=',0,'NaN',1,'Inf',1); md = checkfield(md,'fieldname','amr.deviatoricerror_resolution','numel',[1],'>',0,'<',self.hmax,'NaN',1); md = checkfield(md,'fieldname','amr.deviatoricerror_threshold','numel',[1],'>=',0,'<=',1,'NaN',1); + md = checkfield(md,'fieldname','amr.deviatoricerror_groupthreshold','numel',[1],'>=',0,'<=',1,'NaN',1); md = checkfield(md,'fieldname','amr.deviatoricerror_maximum','numel',[1],'>=',0,'NaN',1,'Inf',1); md = checkfield(md,'fieldname','amr.restart','numel',[1],'>=',0,'<=',1,'NaN',1); @@ -131,7 +137,9 @@ fielddisplay(self,'thicknesserror_resolution',['element length when thickness error estimator is used']); fielddisplay(self,'thicknesserror_threshold',['maximum threshold thickness error permitted']); + fielddisplay(self,'thicknesserror_groupthreshold',['maximum group threshold thickness error permitted']); fielddisplay(self,'thicknesserror_maximum',['maximum thickness error permitted']); fielddisplay(self,'deviatoricerror_resolution',['element length when deviatoric stress error estimator is used']); fielddisplay(self,'deviatoricerror_threshold',['maximum threshold deviatoricstress error permitted']); + fielddisplay(self,'deviatoricerror_groupthreshold',['maximum group threshold deviatoricstress error permitted']); fielddisplay(self,'deviatoricerror_maximum',['maximum deviatoricstress error permitted']); fielddisplay(self,'restart',['indicates if ReMesh() will call before first time step']); @@ -152,7 +160,9 @@ WriteData(fid,prefix,'object',self,'class','amr','fieldname','thicknesserror_resolution','format','Double'); WriteData(fid,prefix,'object',self,'class','amr','fieldname','thicknesserror_threshold','format','Double'); + WriteData(fid,prefix,'object',self,'class','amr','fieldname','thicknesserror_groupthreshold','format','Double'); WriteData(fid,prefix,'object',self,'class','amr','fieldname','thicknesserror_maximum','format','Double'); WriteData(fid,prefix,'object',self,'class','amr','fieldname','deviatoricerror_resolution','format','Double'); WriteData(fid,prefix,'object',self,'class','amr','fieldname','deviatoricerror_threshold','format','Double'); + WriteData(fid,prefix,'object',self,'class','amr','fieldname','deviatoricerror_groupthreshold','format','Double'); WriteData(fid,prefix,'object',self,'class','amr','fieldname','deviatoricerror_maximum','format','Double'); WriteData(fid,prefix,'object',self,'class','amr','fieldname','restart','format','Integer'); Index: /issm/trunk-jpl/src/m/classes/amr.py =================================================================== --- /issm/trunk-jpl/src/m/classes/amr.py (revision 22293) +++ /issm/trunk-jpl/src/m/classes/amr.py (revision 22294) @@ -12,21 +12,23 @@ def __init__(self): # {{{ - self.hmin = 0. - self.hmax = 0. - self.fieldname ='' - self.err = 0. - self.keepmetric = 0. - self.gradation = 0. - self.groundingline_resolution = 0. - self.groundingline_distance = 0. - self.icefront_resolution = 0. - self.icefront_distance = 0. - self.thicknesserror_resolution = 0. - self.thicknesserror_threshold = 0. - self.thicknesserror_maximum = 0. - self.deviatoricerror_resolution= 0. - self.deviatoricerror_threshold = 0. - self.deviatoricerror_maximum = 0. - self.restart=0. + self.hmin = 0. + self.hmax = 0. + self.fieldname ='' + self.err = 0. + self.keepmetric = 0. + self.gradation = 0. + self.groundingline_resolution = 0. + self.groundingline_distance = 0. + self.icefront_resolution = 0. + self.icefront_distance = 0. + self.thicknesserror_resolution = 0. + self.thicknesserror_threshold = 0. + self.thicknesserror_groupthreshold = 0. + self.thicknesserror_maximum = 0. + self.deviatoricerror_resolution = 0. + self.deviatoricerror_threshold = 0. + self.deviatoricerror_groupthreshold = 0. + self.deviatoricerror_maximum = 0. + self.restart = 0. #set defaults self.setdefaultparameters() @@ -45,29 +47,33 @@ string="%s\n%s"%(string,fielddisplay(self,"thicknesserror_resolution","element length when thickness error estimator is used")) string="%s\n%s"%(string,fielddisplay(self,"thicknesserror_threshold","maximum threshold thickness error permitted")) + string="%s\n%s"%(string,fielddisplay(self,"thicknesserror_groupthreshold","maximum group threshold thickness error permitted")) string="%s\n%s"%(string,fielddisplay(self,"thicknesserror_maximum","maximum thickness error permitted")) string="%s\n%s"%(string,fielddisplay(self,"deviatoricerror_resolution","element length when deviatoric stress error estimator is used")) string="%s\n%s"%(string,fielddisplay(self,"deviatoricerror_threshold","maximum threshold deviatoricstress error permitted")) + string="%s\n%s"%(string,fielddisplay(self,"deviatoricerror_groupthreshold","maximum group threshold deviatoric stress error permitted")) string="%s\n%s"%(string,fielddisplay(self,"deviatoricerror_maximum","maximum deviatoricstress error permitted")) - string="%s\n%s"%(string,fielddisplay(self,'restart',['indicates if ReMesh() will call before first time step'])) + string="%s\n%s"%(string,fielddisplay(self,'restart',['indicates if ReMesh() will call before first time step'])) return string #}}} def setdefaultparameters(self): # {{{ - self.hmin = 100. - self.hmax = 100.e3 - self.fieldname = 'Vel' - self.err = 3. - self.keepmetric = 1 - self.gradation = 1.5 - self.groundingline_resolution = 500. - self.groundingline_distance = 0 - self.icefront_resolution = 500. - self.icefront_distance = 0 - self.thicknesserror_resolution = 500. - self.thicknesserror_threshold = 0 - self.thicknesserror_maximum = 0 - self.deviatoricerror_resolution= 500. - self.deviatoricerror_threshold = 0 - self.deviatoricerror_maximum = 0 - self.restart = 0. + self.hmin = 100. + self.hmax = 100.e3 + self.fieldname = 'Vel' + self.err = 3. + self.keepmetric = 1 + self.gradation = 1.5 + self.groundingline_resolution = 500. + self.groundingline_distance = 0 + self.icefront_resolution = 500. + self.icefront_distance = 0 + self.thicknesserror_resolution = 500. + self.thicknesserror_threshold = 0 + self.thicknesserror_groupthreshold = 0 + self.thicknesserror_maximum = 0 + self.deviatoricerror_resolution = 500. + self.deviatoricerror_threshold = 0 + self.deviatoricerror_groupthreshold = 0 + self.deviatoricerror_maximum = 0 + self.restart = 0. return self #}}} @@ -83,9 +89,11 @@ md = checkfield(md,'fieldname','amr.thicknesserror_resolution','numel',[1],'>',0,'<',self.hmax,'NaN',1); md = checkfield(md,'fieldname','amr.thicknesserror_threshold','numel',[1],'>=',0,'<=',1,'NaN',1); + md = checkfield(md,'fieldname','amr.thicknesserror_groupthreshold','numel',[1],'>=',0,'<=',1,'NaN',1); md = checkfield(md,'fieldname','amr.thicknesserror_maximum','numel',[1],'>=',0,'NaN',1,'Inf',1); md = checkfield(md,'fieldname','amr.deviatoricerror_resolution','numel',[1],'>',0,'<',self.hmax,'NaN',1); md = checkfield(md,'fieldname','amr.deviatoricerror_threshold','numel',[1],'>=',0,'<=',1,'NaN',1); + md = checkfield(md,'fieldname','amr.deviatoricerror_groupthreshold','numel',[1],'>=',0,'<=',1,'NaN',1); md = checkfield(md,'fieldname','amr.deviatoricerror_maximum','numel',[1],'>=',0,'NaN',1,'Inf',1); - md = checkfield(md,'fieldname','amr.restart','numel',[1],'>=',0,'<=',1,'NaN',1) + md = checkfield(md,'fieldname','amr.restart','numel',[1],'>=',0,'<=',1,'NaN',1) return md # }}} @@ -104,8 +112,10 @@ WriteData(fid,prefix,'object',self,'fieldname','thicknesserror_resolution','format','Double'); WriteData(fid,prefix,'object',self,'fieldname','thicknesserror_threshold','format','Double'); + WriteData(fid,prefix,'object',self,'fieldname','thicknesserror_groupthreshold','format','Double'); WriteData(fid,prefix,'object',self,'fieldname','thicknesserror_maximum','format','Double'); WriteData(fid,prefix,'object',self,'fieldname','deviatoricerror_resolution','format','Double'); WriteData(fid,prefix,'object',self,'fieldname','deviatoricerror_threshold','format','Double'); + WriteData(fid,prefix,'object',self,'fieldname','deviatoricerror_groupthreshold','format','Double'); WriteData(fid,prefix,'object',self,'fieldname','deviatoricerror_maximum','format','Double'); - WriteData(fid,prefix,'object',self,'class','amr','fieldname','restart','format','Integer') + WriteData(fid,prefix,'object',self,'class','amr','fieldname','restart','format','Integer') # }}} Index: /issm/trunk-jpl/src/m/classes/amrneopz.m =================================================================== --- /issm/trunk-jpl/src/m/classes/amrneopz.m (revision 22293) +++ /issm/trunk-jpl/src/m/classes/amrneopz.m (revision 22294) @@ -6,14 +6,16 @@ classdef amrneopz properties (SetAccess=public) - level_max = 0; - gradation = 0; - lag = 0; - groundingline_distance = 0; - icefront_distance = 0; - thicknesserror_threshold = 0; - thicknesserror_maximum = 0; - deviatoricerror_threshold = 0; - deviatoricerror_maximum = 0; - restart = 0; + level_max = 0; + gradation = 0; + lag = 0; + groundingline_distance = 0; + icefront_distance = 0; + thicknesserror_threshold = 0; + thicknesserror_groupthreshold = 0; + thicknesserror_maximum = 0; + deviatoricerror_threshold = 0; + deviatoricerror_groupthreshold = 0; + deviatoricerror_maximum = 0; + restart = 0; end methods (Static) @@ -59,4 +61,7 @@ %radius_level_max was deleted! + %2017 December 18th + %group threshold was inserted + %end% }}} end @@ -73,16 +78,18 @@ %filter parameters: - self.level_max = 2; - self.gradation = 1.5; - self.lag = 1.1; + self.level_max = 2; + self.gradation = 1.5; + self.lag = 1.1; %other criterias - self.groundingline_distance = 10000; - self.icefront_distance = 0; - self.thicknesserror_threshold = 0; - self.thicknesserror_maximum = 0; - self.deviatoricerror_threshold= 0; - self.deviatoricerror_maximum = 0; - self.restart = 0; + self.groundingline_distance = 10000; + self.icefront_distance = 0; + self.thicknesserror_threshold = 0; + self.thicknesserror_groupthreshold = 0; + self.thicknesserror_maximum = 0; + self.deviatoricerror_threshold = 0; + self.deviatoricerror_groupthreshold = 0; + self.deviatoricerror_maximum = 0; + self.restart = 0; end % }}} function md = checkconsistency(self,md,solution,analyses) % {{{ @@ -94,6 +101,8 @@ md = checkfield(md,'fieldname','amr.icefront_distance','numel',[1],'>=',0,'NaN',1,'Inf',1); md = checkfield(md,'fieldname','amr.thicknesserror_threshold','numel',[1],'>=',0,'<=',1,'NaN',1); + md = checkfield(md,'fieldname','amr.thicknesserror_groupthreshold','numel',[1],'>=',0,'<=',1,'NaN',1); md = checkfield(md,'fieldname','amr.thicknesserror_maximum','numel',[1],'>=',0,'NaN',1,'Inf',1); md = checkfield(md,'fieldname','amr.deviatoricerror_threshold','numel',[1],'>=',0,'<=',1,'NaN',1); + md = checkfield(md,'fieldname','amr.deviatoricerror_groupthreshold','numel',[1],'>=',0,'<=',1,'NaN',1); md = checkfield(md,'fieldname','amr.deviatoricerror_maximum','numel',[1],'>=',0,'NaN',1,'Inf',1); md = checkfield(md,'fieldname','amr.restart','numel',[1],'>=',0,'<=',1,'NaN',1); @@ -108,6 +117,8 @@ fielddisplay(self,'icefront_distance',['distance around the ice front which elements will be refined']); fielddisplay(self,'thicknesserror_threshold',['maximum threshold thickness error permitted']); + fielddisplay(self,'thicknesserror_groupthreshold',['maximum group threshold thickness error permitted']); fielddisplay(self,'thicknesserror_maximum',['maximum thickness error permitted']); fielddisplay(self,'deviatoricerror_threshold',['maximum threshold deviatoricstress error permitted']); + fielddisplay(self,'deviatoricerror_groupthreshold',['maximum group threshold deviatoricstress error permitted']); fielddisplay(self,'deviatoricerror_maximum',['maximum deviatoricstress error permitted']); fielddisplay(self,'restart',['indicates if ReMesh() will call before first time step']); @@ -122,6 +133,8 @@ WriteData(fid,prefix,'object',self,'class','amr','fieldname','icefront_distance','format','Double'); WriteData(fid,prefix,'object',self,'class','amr','fieldname','thicknesserror_threshold','format','Double'); + WriteData(fid,prefix,'object',self,'class','amr','fieldname','thicknesserror_groupthreshold','format','Double'); WriteData(fid,prefix,'object',self,'class','amr','fieldname','thicknesserror_maximum','format','Double'); WriteData(fid,prefix,'object',self,'class','amr','fieldname','deviatoricerror_threshold','format','Double'); + WriteData(fid,prefix,'object',self,'class','amr','fieldname','deviatoricerror_groupthreshold','format','Double'); WriteData(fid,prefix,'object',self,'class','amr','fieldname','deviatoricerror_maximum','format','Double'); WriteData(fid,prefix,'object',self,'class','amr','fieldname','restart','format','Integer'); @@ -135,6 +148,8 @@ writejsdouble(fid,[modelname '.amr.icefront_distance'],self.icefront_distance); writejsdouble(fid,[modelname '.amr.thicknesserror_threshold'],self.thicknesserror_threshold); + writejsdouble(fid,[modelname '.amr.thicknesserror_groupthreshold'],self.thicknesserror_threshold); writejsdouble(fid,[modelname '.amr.thicknesserror_maximum'],self.thicknesserror_maximum); writejsdouble(fid,[modelname '.amr.deviatoricerror_threshold'],self.deviatoricerror_threshold); + writejsdouble(fid,[modelname '.amr.deviatoricerror_groupthreshold'],self.deviatoricerror_threshold); writejsdouble(fid,[modelname '.amr.deviatoricerror_maximum'],self.deviatoricerror_maximum); writejsdouble(fid,[modelname '.amr.restart'],self.restart); Index: /issm/trunk-jpl/src/m/contrib/tsantos/mismip/ice_evolution.m =================================================================== --- /issm/trunk-jpl/src/m/contrib/tsantos/mismip/ice_evolution.m (revision 22293) +++ /issm/trunk-jpl/src/m/contrib/tsantos/mismip/ice_evolution.m (revision 22294) @@ -7,15 +7,15 @@ % % Default: y=40km, i0=1 -% [ga iv ivaf GL_y nelem t] = ice_evolution(md); +% [ga iv ivaf GL_y IF_y nelem t] = ice_evolution(md); % % Default: y=40km -% [ga iv ivaf GL_y nelem t] = ice_evolution(md,i0); +% [ga iv ivaf GL_y IF_y nelem t] = ice_evolution(md,i0); % % Use this for y the borders (y=0 or y=ymax) -% [ga iv ivaf GL_y nelem t] = ice_evolution(md,i0,y); +% [ga iv ivaf GL_y IF_y nelem t] = ice_evolution(md,i0,y); % % -function [ga iv ivaf GL_y nelem t] = ice_evolution(varargin), +function [ga iv ivaf GL_y IF_y nelem t] = ice_evolution(varargin), ga = []; @@ -23,4 +23,5 @@ ivaf = []; GL_y = []; + IF_y = []; nelem = []; t = []; @@ -66,5 +67,5 @@ end t(i) = md.results.TransientSolution(i).time; - %find GL position between y0 and y1 + %find GL position between y0 and y1 {{{ [glx gly] = gl_position(md,i,0); pos = find(glyy0); @@ -87,5 +88,30 @@ end end - + %}}} + %find IF position between y0 and y1 {{{ + [ifx ify] = if_position(md,i,0); + if(length(ifx)==0) + continue + end + pos = find(ifyy0); + x = ify(pos); + v = ifx(pos); + if(length(pos)==0) + error('pos is null') + elseif(length(pos)==1) + %this should be used for y=0 or y=ymax + IF_y(i) = v; + else + %this should be used when y is inside the domain; so, use linear interpolation + xq = [y0:dy:y1]; + vq = interp1(x,v,xq,'linear'); + pos = find(xq==y); + if(pos) + IF_y(i) = vq(pos); + else + error('pos is null') + end + end + %}}} end