Changeset 24048


Ignore:
Timestamp:
06/25/19 15:00:14 (6 years ago)
Author:
youngmc3
Message:

CHG: added CalvingFluxLevelset

Location:
issm/trunk-jpl/src/c
Files:
10 edited

Legend:

Unmodified
Added
Removed

  • issm/trunk-jpl/src/c/classes/Elements/Element.cpp

    r24045 r24048  
    32923292                        }
    32933293                        break;
     3294                case CalvingFluxLevelsetEnum: this->CalvingFluxLevelset(); break;
    32943295                case StrainRateparallelEnum: this->StrainRateparallel(); break;
    32953296                case StrainRateperpendicularEnum: this->StrainRateperpendicular(); break;

  • issm/trunk-jpl/src/c/classes/Elements/Element.h

    r24020 r24048  
    208208                virtual void       CalvingCrevasseDepth(void){_error_("not implemented yet");};
    209209                virtual void        CalvingRateLevermann(void)=0;
     210                virtual void       CalvingFluxLevelset(void){_error_("not implemented yet");};
    210211                virtual IssmDouble CharacteristicLength(void)=0;
    211212                virtual void       ComputeBasalStress(void){_error_("not implemented yet");};

  • issm/trunk-jpl/src/c/classes/Elements/Tria.cpp

    r24020 r24048  
    482482        delete gauss;
    483483
     484}
     485/*}}}*/
     486void       Tria::CalvingFluxLevelset(){/*{{{*/
     487
     488        /*Make sure there is an ice front here*/
     489        if(!IsIceInElement() || !IsZeroLevelset(MaskIceLevelsetEnum)){
     490                IssmDouble flux_per_area=0;
     491                this->inputs->AddInput(new TriaInput(CalvingFluxLevelsetEnum,&flux_per_area,P0Enum));
     492        }
     493        else{
     494        int               domaintype,index1,index2;
     495        const IssmPDouble epsilon = 1.e-15;
     496        IssmDouble        s1,s2;
     497        IssmDouble        gl[NUMVERTICES];
     498        IssmDouble        xyz_front[2][3];
     499
     500
     501        IssmDouble *xyz_list = NULL;
     502        this->GetVerticesCoordinates(&xyz_list);
     503
     504        /*Recover parameters and values*/
     505        parameters->FindParam(&domaintype,DomainTypeEnum);
     506        GetInputListOnVertices(&gl[0],MaskIceLevelsetEnum);
     507
     508        /*Be sure that values are not zero*/
     509        if(gl[0]==0.) gl[0]=gl[0]+epsilon;
     510        if(gl[1]==0.) gl[1]=gl[1]+epsilon;
     511        if(gl[2]==0.) gl[2]=gl[2]+epsilon;
     512
     513        if(domaintype==Domain2DverticalEnum){
     514                _error_("not implemented");
     515        }
     516        else if(domaintype==Domain2DhorizontalEnum || domaintype==Domain3DEnum || domaintype==Domain3DsurfaceEnum){
     517                int pt1 = 0;
     518                int pt2 = 1;
     519                if(gl[0]*gl[1]>0){ //Nodes 0 and 1 are similar, so points must be found on segment 0-2 and 1-2
     520
     521                        /*Portion of the segments*/
     522                        s1=gl[2]/(gl[2]-gl[1]);
     523                        s2=gl[2]/(gl[2]-gl[0]);
     524                        if(gl[2]<0.){
     525                                pt1 = 1; pt2 = 0;
     526                        }
     527                        xyz_front[pt2][0]=xyz_list[3*2+0]+s1*(xyz_list[3*1+0]-xyz_list[3*2+0]);
     528                        xyz_front[pt2][1]=xyz_list[3*2+1]+s1*(xyz_list[3*1+1]-xyz_list[3*2+1]);
     529                        xyz_front[pt2][2]=xyz_list[3*2+2]+s1*(xyz_list[3*1+2]-xyz_list[3*2+2]);
     530                        xyz_front[pt1][0]=xyz_list[3*2+0]+s2*(xyz_list[3*0+0]-xyz_list[3*2+0]);
     531                        xyz_front[pt1][1]=xyz_list[3*2+1]+s2*(xyz_list[3*0+1]-xyz_list[3*2+1]);
     532                        xyz_front[pt1][2]=xyz_list[3*2+2]+s2*(xyz_list[3*0+2]-xyz_list[3*2+2]);
     533                }
     534                else if(gl[1]*gl[2]>0){ //Nodes 1 and 2 are similar, so points must be found on segment 0-1 and 0-2
     535
     536                        /*Portion of the segments*/
     537                        s1=gl[0]/(gl[0]-gl[1]);
     538                        s2=gl[0]/(gl[0]-gl[2]);
     539                        if(gl[0]<0.){
     540                                pt1 = 1; pt2 = 0;
     541                        }
     542
     543                        xyz_front[pt1][0]=xyz_list[3*0+0]+s1*(xyz_list[3*1+0]-xyz_list[3*0+0]);
     544                        xyz_front[pt1][1]=xyz_list[3*0+1]+s1*(xyz_list[3*1+1]-xyz_list[3*0+1]);
     545                        xyz_front[pt1][2]=xyz_list[3*0+2]+s1*(xyz_list[3*1+2]-xyz_list[3*0+2]);
     546                        xyz_front[pt2][0]=xyz_list[3*0+0]+s2*(xyz_list[3*2+0]-xyz_list[3*0+0]);
     547                        xyz_front[pt2][1]=xyz_list[3*0+1]+s2*(xyz_list[3*2+1]-xyz_list[3*0+1]);
     548                        xyz_front[pt2][2]=xyz_list[3*0+2]+s2*(xyz_list[3*2+2]-xyz_list[3*0+2]);
     549                }
     550                else if(gl[0]*gl[2]>0){ //Nodes 0 and 2 are similar, so points must be found on segment 1-0 and 1-2
     551
     552                        /*Portion of the segments*/
     553                        s1=gl[1]/(gl[1]-gl[0]);
     554                        s2=gl[1]/(gl[1]-gl[2]);
     555                        if(gl[1]<0.){
     556                                pt1 = 1; pt2 = 0;
     557                        }
     558
     559                        xyz_front[pt2][0]=xyz_list[3*1+0]+s1*(xyz_list[3*0+0]-xyz_list[3*1+0]);
     560                        xyz_front[pt2][1]=xyz_list[3*1+1]+s1*(xyz_list[3*0+1]-xyz_list[3*1+1]);
     561                        xyz_front[pt2][2]=xyz_list[3*1+2]+s1*(xyz_list[3*0+2]-xyz_list[3*1+2]);
     562                        xyz_front[pt1][0]=xyz_list[3*1+0]+s2*(xyz_list[3*2+0]-xyz_list[3*1+0]);
     563                        xyz_front[pt1][1]=xyz_list[3*1+1]+s2*(xyz_list[3*2+1]-xyz_list[3*1+1]);
     564                        xyz_front[pt1][2]=xyz_list[3*1+2]+s2*(xyz_list[3*2+2]-xyz_list[3*1+2]);
     565                }
     566                else{
     567                        _error_("case not possible");
     568                }
     569
     570        }
     571        else _error_("mesh type "<<EnumToStringx(domaintype)<<"not supported yet ");
     572
     573        /*Some checks in debugging mode*/
     574        _assert_(s1>=0 && s1<=1.);
     575        _assert_(s2>=0 && s2<=1.);
     576
     577        /*Get normal vector*/
     578        IssmDouble normal[3];
     579        this->NormalSection(&normal[0],&xyz_front[0][0]);
     580        normal[0] = -normal[0];
     581        normal[1] = -normal[1];
     582
     583        /*Get inputs*/
     584        IssmDouble flux = 0.;
     585        IssmDouble area = 0.;
     586        IssmDouble calvingratex,calvingratey,thickness,Jdet,flux_per_area;
     587        IssmDouble rho_ice=FindParam(MaterialsRhoIceEnum);
     588        Input* thickness_input=inputs->GetInput(ThicknessEnum); _assert_(thickness_input);
     589        Input* calvingratex_input=NULL;
     590        Input* calvingratey_input=NULL;
     591        if(domaintype==Domain2DhorizontalEnum){
     592                calvingratex_input=inputs->GetInput(CalvingratexEnum); _assert_(calvingratex_input);
     593                calvingratey_input=inputs->GetInput(CalvingrateyEnum); _assert_(calvingratey_input);
     594        }
     595        else{
     596                calvingratex_input=inputs->GetInput(CalvingratexAverageEnum); _assert_(calvingratex_input);
     597                calvingratey_input=inputs->GetInput(CalvingrateyAverageEnum); _assert_(calvingratey_input);
     598        }
     599
     600        /*Start looping on Gaussian points*/
     601        Gauss* gauss=this->NewGauss(xyz_list,&xyz_front[0][0],3);
     602        for(int ig=gauss->begin();ig<gauss->end();ig++){
     603
     604                gauss->GaussPoint(ig);
     605                thickness_input->GetInputValue(&thickness,gauss);
     606                calvingratex_input->GetInputValue(&calvingratex,gauss);
     607                calvingratey_input->GetInputValue(&calvingratey,gauss);
     608                this->JacobianDeterminantSurface(&Jdet,&xyz_front[0][0],gauss);
     609
     610                flux += rho_ice*Jdet*gauss->weight*thickness*(calvingratex*normal[0] + calvingratey*normal[1]);
     611                area += Jdet*gauss->weight*thickness;
     612
     613                flux_per_area=flux/area;
     614        }
     615       
     616        this->inputs->AddInput(new TriaInput(CalvingFluxLevelsetEnum,&flux_per_area,P0Enum));
     617        }
    484618}
    485619/*}}}*/

  • issm/trunk-jpl/src/c/classes/Elements/Tria.h

    r24020 r24048  
    5454                void        CalvingCrevasseDepth();
    5555                void                    CalvingRateLevermann();
     56                void                    CalvingFluxLevelset();
    5657                IssmDouble  CharacteristicLength(void);
    5758                void        ComputeBasalStress(void);

  • issm/trunk-jpl/src/c/classes/FemModel.cpp

    r24047 r24048  
    10291029                Element* element=dynamic_cast<Element*>(this->elements->GetObjectByOffset(i));
    10301030                element->CalvingRateLevermann();
     1031        }
     1032}
     1033/*}}}*/
     1034void FemModel::CalvingFluxLevelsetx(){/*{{{*/
     1035
     1036        for(int i=0;i<elements->Size();i++){
     1037                Element* element=dynamic_cast<Element*>(this->elements->GetObjectByOffset(i));
     1038                element->CalvingFluxLevelset();
    10311039        }
    10321040}

  • issm/trunk-jpl/src/c/classes/FemModel.h

    r24047 r24048  
    9090                void CalvingRateVonmisesx();
    9191                void CalvingRateLevermannx();
     92                void CalvingFluxLevelsetx();
    9293                void DeviatoricStressx();
    9394                void Divergencex(IssmDouble* pdiv);

  • issm/trunk-jpl/src/c/shared/Enum/Enum.vim

    r24045 r24048  
    482482syn keyword cConstant CalvingrateyAverageEnum
    483483syn keyword cConstant CalvingrateyEnum
     484syn keyword cConstant CalvingFluxLevelsetEnum
    484485syn keyword cConstant ConvergedEnum
    485486syn keyword cConstant CrevasseDepthEnum
    … …  
    12711272syn keyword cType Cfsurfacesquare
    12721273syn keyword cType Channel
     1274syn keyword cType classes
    12731275syn keyword cType Constraint
    12741276syn keyword cType Constraints
    … …  
    12771279syn keyword cType ControlInput
    12781280syn keyword cType Covertree
     1281syn keyword cType DatasetInput
    12791282syn keyword cType DataSetParam
    1280 syn keyword cType DatasetInput
    12811283syn keyword cType Definition
    12821284syn keyword cType DependentObject
    … …  
    12911293syn keyword cType ElementHook
    12921294syn keyword cType ElementMatrix
     1295syn keyword cType Elements
    12931296syn keyword cType ElementVector
    1294 syn keyword cType Elements
    12951297syn keyword cType ExponentialVariogram
    12961298syn keyword cType ExternalResult
    … …  
    12991301syn keyword cType Friction
    13001302syn keyword cType Gauss
     1303syn keyword cType GaussianVariogram
     1304syn keyword cType gaussobjects
    13011305syn keyword cType GaussPenta
    13021306syn keyword cType GaussSeg
    13031307syn keyword cType GaussTetra
    13041308syn keyword cType GaussTria
    1305 syn keyword cType GaussianVariogram
    13061309syn keyword cType GenericExternalResult
    13071310syn keyword cType GenericOption
    … …  
    13181321syn keyword cType IssmDirectApplicInterface
    13191322syn keyword cType IssmParallelDirectApplicInterface
     1323syn keyword cType krigingobjects
    13201324syn keyword cType Load
    13211325syn keyword cType Loads
    … …  
    13281332syn keyword cType Matice
    13291333syn keyword cType Matlitho
     1334syn keyword cType matrixobjects
    13301335syn keyword cType MatrixParam
    13311336syn keyword cType Misfit
    … …  
    13401345syn keyword cType Observations
    13411346syn keyword cType Option
     1347syn keyword cType Options
    13421348syn keyword cType OptionUtilities
    1343 syn keyword cType Options
    13441349syn keyword cType Param
    13451350syn keyword cType Parameters
    … …  
    13541359syn keyword cType Regionaloutput
    13551360syn keyword cType Results
     1361syn keyword cType Riftfront
    13561362syn keyword cType RiftStruct
    1357 syn keyword cType Riftfront
    13581363syn keyword cType Seg
    13591364syn keyword cType SegInput
     1365syn keyword cType Segment
    13601366syn keyword cType SegRef
    1361 syn keyword cType Segment
    13621367syn keyword cType SpcDynamic
    13631368syn keyword cType SpcStatic
    … …  
    13791384syn keyword cType Vertex
    13801385syn keyword cType Vertices
    1381 syn keyword cType classes
    1382 syn keyword cType gaussobjects
    1383 syn keyword cType krigingobjects
    1384 syn keyword cType matrixobjects
    13851386syn keyword cType AdjointBalancethickness2Analysis
    13861387syn keyword cType AdjointBalancethicknessAnalysis
    … …  
    13911392syn keyword cType BalancethicknessSoftAnalysis
    13921393syn keyword cType BalancevelocityAnalysis
     1394syn keyword cType DamageCalvingAnalysis
    13931395syn keyword cType DamageEvolutionAnalysis
    13941396syn keyword cType DepthAverageAnalysis
    … …  
    14011403syn keyword cType FreeSurfaceBaseAnalysis
    14021404syn keyword cType FreeSurfaceTopAnalysis
     1405syn keyword cType GiaIvinsAnalysis
    14031406syn keyword cType GLheightadvectionAnalysis
    1404 syn keyword cType GiaIvinsAnalysis
    14051407syn keyword cType HydrologyDCEfficientAnalysis
    14061408syn keyword cType HydrologyDCInefficientAnalysis

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

    r24045 r24048  
    478478        CalvingrateyAverageEnum,
    479479        CalvingrateyEnum,
     480        CalvingFluxLevelsetEnum,
    480481        ConvergedEnum,
    481482        CrevasseDepthEnum,

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

    r24045 r24048  
    484484                case CalvingrateyAverageEnum : return "CalvingrateyAverage";
    485485                case CalvingrateyEnum : return "Calvingratey";
     486                case CalvingFluxLevelsetEnum : return "CalvingFluxLevelset";
    486487                case ConvergedEnum : return "Converged";
    487488                case CrevasseDepthEnum : return "CrevasseDepth";

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

    r24045 r24048  
    493493              else if (strcmp(name,"CalvingrateyAverage")==0) return CalvingrateyAverageEnum;
    494494              else if (strcmp(name,"Calvingratey")==0) return CalvingrateyEnum;
     495              else if (strcmp(name,"CalvingFluxLevelset")==0) return CalvingFluxLevelsetEnum;
    495496              else if (strcmp(name,"Converged")==0) return ConvergedEnum;
    496497              else if (strcmp(name,"CrevasseDepth")==0) return CrevasseDepthEnum;
    … …  
    505506              else if (strcmp(name,"DeviatoricStressyy")==0) return DeviatoricStressyyEnum;
    506507              else if (strcmp(name,"DeviatoricStressyz")==0) return DeviatoricStressyzEnum;
    507               else if (strcmp(name,"DeviatoricStresszz")==0) return DeviatoricStresszzEnum;
    508508         else stage=5;
    509509   }
    510510   if(stage==5){
    511               if (strcmp(name,"DeviatoricStress1")==0) return DeviatoricStress1Enum;
     511              if (strcmp(name,"DeviatoricStresszz")==0) return DeviatoricStresszzEnum;
     512              else if (strcmp(name,"DeviatoricStress1")==0) return DeviatoricStress1Enum;
    512513              else if (strcmp(name,"DeviatoricStress2")==0) return DeviatoricStress2Enum;
    513514              else if (strcmp(name,"DistanceToCalvingfront")==0) return DistanceToCalvingfrontEnum;
    … …  
    628629              else if (strcmp(name,"Misfit")==0) return MisfitEnum;
    629630              else if (strcmp(name,"Neumannflux")==0) return NeumannfluxEnum;
    630               else if (strcmp(name,"NewDamage")==0) return NewDamageEnum;
    631631         else stage=6;
    632632   }
    633633   if(stage==6){
    634               if (strcmp(name,"Node")==0) return NodeEnum;
     634              if (strcmp(name,"NewDamage")==0) return NewDamageEnum;
     635              else if (strcmp(name,"Node")==0) return NodeEnum;
    635636              else if (strcmp(name,"OmegaAbsGradient")==0) return OmegaAbsGradientEnum;
    636637              else if (strcmp(name,"P0")==0) return P0Enum;
    … …  
    751752              else if (strcmp(name,"SmbVz")==0) return SmbVzEnum;
    752753              else if (strcmp(name,"SmbW")==0) return SmbWEnum;
    753               else if (strcmp(name,"SmbWini")==0) return SmbWiniEnum;
    754754         else stage=7;
    755755   }
    756756   if(stage==7){
    757               if (strcmp(name,"SmbZMax")==0) return SmbZMaxEnum;
     757              if (strcmp(name,"SmbWini")==0) return SmbWiniEnum;
     758              else if (strcmp(name,"SmbZMax")==0) return SmbZMaxEnum;
    758759              else if (strcmp(name,"SmbZMin")==0) return SmbZMinEnum;
    759760              else if (strcmp(name,"SmbZTop")==0) return SmbZTopEnum;
    … …  
    874875              else if (strcmp(name,"Outputdefinition59")==0) return Outputdefinition59Enum;
    875876              else if (strcmp(name,"Outputdefinition5")==0) return Outputdefinition5Enum;
    876               else if (strcmp(name,"Outputdefinition60")==0) return Outputdefinition60Enum;
    877877         else stage=8;
    878878   }
    879879   if(stage==8){
    880               if (strcmp(name,"Outputdefinition61")==0) return Outputdefinition61Enum;
     880              if (strcmp(name,"Outputdefinition60")==0) return Outputdefinition60Enum;
     881              else if (strcmp(name,"Outputdefinition61")==0) return Outputdefinition61Enum;
    881882              else if (strcmp(name,"Outputdefinition62")==0) return Outputdefinition62Enum;
    882883              else if (strcmp(name,"Outputdefinition63")==0) return Outputdefinition63Enum;
    … …  
    997998              else if (strcmp(name,"DoubleVecParam")==0) return DoubleVecParamEnum;
    998999              else if (strcmp(name,"Element")==0) return ElementEnum;
    999               else if (strcmp(name,"ElementHook")==0) return ElementHookEnum;
    10001000         else stage=9;
    10011001   }
    10021002   if(stage==9){
    1003               if (strcmp(name,"ElementSId")==0) return ElementSIdEnum;
     1003              if (strcmp(name,"ElementHook")==0) return ElementHookEnum;
     1004              else if (strcmp(name,"ElementSId")==0) return ElementSIdEnum;
    10041005              else if (strcmp(name,"EnthalpyAnalysis")==0) return EnthalpyAnalysisEnum;
    10051006              else if (strcmp(name,"EsaAnalysis")==0) return EsaAnalysisEnum;
    … …  
    11201121              else if (strcmp(name,"MaxAbsVy")==0) return MaxAbsVyEnum;
    11211122              else if (strcmp(name,"MaxAbsVz")==0) return MaxAbsVzEnum;
    1122               else if (strcmp(name,"MaxDivergence")==0) return MaxDivergenceEnum;
    11231123         else stage=10;
    11241124   }
    11251125   if(stage==10){
    1126               if (strcmp(name,"MaxVel")==0) return MaxVelEnum;
     1126              if (strcmp(name,"MaxDivergence")==0) return MaxDivergenceEnum;
     1127              else if (strcmp(name,"MaxVel")==0) return MaxVelEnum;
    11271128              else if (strcmp(name,"MaxVx")==0) return MaxVxEnum;
    11281129              else if (strcmp(name,"MaxVy")==0) return MaxVyEnum;
    … …  
    12431244              else if (strcmp(name,"Tetra")==0) return TetraEnum;
    12441245              else if (strcmp(name,"TetraInput")==0) return TetraInputEnum;
    1245               else if (strcmp(name,"ThermalAnalysis")==0) return ThermalAnalysisEnum;
    12461246         else stage=11;
    12471247   }
    12481248   if(stage==11){
    1249               if (strcmp(name,"ThermalSolution")==0) return ThermalSolutionEnum;
     1249              if (strcmp(name,"ThermalAnalysis")==0) return ThermalAnalysisEnum;
     1250              else if (strcmp(name,"ThermalSolution")==0) return ThermalSolutionEnum;
    12501251              else if (strcmp(name,"ThicknessErrorEstimator")==0) return ThicknessErrorEstimatorEnum;
    12511252              else if (strcmp(name,"TotalCalvingFluxLevelset")==0) return TotalCalvingFluxLevelsetEnum;
Note: See TracChangeset for help on using the changeset viewer.