Changeset 22481


Ignore:
Timestamp:
02/27/18 14:28:51 (7 years ago)
Author:
Eric.Larour
Message:

CHG: implemented basal forcings in vectorial format.

Location:
issm/branches/trunk-larour-NatGeoScience2016/src
Files:
1 added
26 edited

Legend:

Unmodified
Added
Removed

  • issm/branches/trunk-larour-NatGeoScience2016/src/c/analyses/MasstransportAnalysis.cpp

    r22143 r22481  
    123123        bool   ismovingfront;
    124124        bool   issmb;
     125        int    basalforcingsmodel;
    125126
    126127        /*Fetch data needed: */
    … …  
    130131        iomodel->FindConstant(&ismovingfront,"md.transient.ismovingfront");
    131132        iomodel->FindConstant(&issmb,"md.transient.issmb");
     133        iomodel->FindConstant(&basalforcingsmodel,"md.basalforcings.model");
    132134
    133135        /*Finite element type*/
    … …  
    155157        iomodel->FetchDataToInput(elements,"md.basalforcings.groundedice_melting_rate",BasalforcingsGroundediceMeltingRateEnum);
    156158        iomodel->FetchDataToInput(elements,"md.basalforcings.floatingice_melting_rate",BasalforcingsFloatingiceMeltingRateEnum);
     159        if (basalforcingsmodel==SpatialLinearFloatingMeltRateEnum){
     160                iomodel->FetchDataToInput(elements,"md.basalforcings.deepwater_melting_rate",BasalforcingsDeepwaterMeltingRateEnum);
     161                iomodel->FetchDataToInput(elements,"md.basalforcings.deepwater_elevation",BasalforcingsDeepwaterElevationEnum);
     162                iomodel->FetchDataToInput(elements,"md.basalforcings.upperwater_elevation",BasalforcingsUpperwaterElevationEnum);
     163        }
    157164        iomodel->FetchDataToInput(elements,"md.initialization.vx",VxEnum);
    158165        iomodel->FetchDataToInput(elements,"md.initialization.vy",VyEnum);

  • issm/branches/trunk-larour-NatGeoScience2016/src/c/analyses/ThermalAnalysis.cpp

    r21759 r22481  
    754754        IssmDouble *xyz_list  = NULL;
    755755        IssmDouble  n=3.0;
    756         bool        hack      = false;
     756        bool        hack      = true;
    757757
    758758        /*Fetch number of nodes and dof for this finite element*/

  • issm/branches/trunk-larour-NatGeoScience2016/src/c/classes/Elements/Element.cpp

    r22248 r22481  
    17301730
    17311731}/*}}}*/
     1732void       Element::SpatialLinearFloatingiceMeltingRate(){/*{{{*/
     1733
     1734        int numvertices      = this->GetNumberOfVertices();
     1735        IssmDouble* deepwatermelt     = xNew<IssmDouble>(numvertices);
     1736        IssmDouble* deepwaterel     = xNew<IssmDouble>(numvertices);
     1737        IssmDouble* upperwaterel     = xNew<IssmDouble>(numvertices);
     1738        IssmDouble* base     = xNew<IssmDouble>(numvertices);
     1739        IssmDouble* values   = xNew<IssmDouble>(numvertices);
     1740
     1741        this->GetInputListOnVertices(base,BaseEnum);
     1742        this->GetInputListOnVertices(deepwatermelt,BasalforcingsDeepwaterMeltingRateEnum);
     1743        this->GetInputListOnVertices(deepwaterel,BasalforcingsDeepwaterElevationEnum);
     1744        this->GetInputListOnVertices(upperwaterel,BasalforcingsUpperwaterElevationEnum);
     1745       
     1746        for(int i=0;i<numvertices;i++){
     1747                if(base[i]>upperwaterel[i])      values[i]=0;
     1748                else if (base[i]<deepwaterel[i]) values[i]=deepwatermelt[i];
     1749                else values[i]=deepwatermelt[i]*(base[i]-upperwaterel[i])/(deepwaterel[i]-upperwaterel[i]);
     1750        }
     1751
     1752        this->AddInput(BasalforcingsFloatingiceMeltingRateEnum,values,P1Enum);
     1753        xDelete<IssmDouble>(base);
     1754        xDelete<IssmDouble>(deepwaterel);
     1755        xDelete<IssmDouble>(deepwatermelt);
     1756        xDelete<IssmDouble>(upperwaterel);
     1757        xDelete<IssmDouble>(values);
     1758
     1759}/*}}}*/
    17321760void       Element::MantlePlumeGeothermalFlux(){/*{{{*/
    17331761

  • issm/branches/trunk-larour-NatGeoScience2016/src/c/classes/Elements/Element.h

    r22114 r22481  
    126126                bool               IsWaterInElement();
    127127                void               LinearFloatingiceMeltingRate();
     128                void               SpatialLinearFloatingiceMeltingRate();
    128129                void               MantlePlumeGeothermalFlux();
    129130                void               MarshallElement(char** pmarshalled_data,int* pmarshalled_data_size, int marshall_direction,int numanalyses);

  • issm/branches/trunk-larour-NatGeoScience2016/src/c/classes/Elements/Tria.cpp

    r22248 r22481  
    43894389        if(!(this->inputs->Max(MaskIceLevelsetEnum)<0) && !IsWaterInElement()) return;
    43904390
     4391        /*early return if we are fully floating: */
     4392        if (this->inputs->Max(MaskGroundediceLevelsetEnum)<=0)return;
     4393
    43914394        /*recover computational flags: */
    43924395        this->parameters->FindParam(&computerigid,SealevelriseRigidEnum);

  • issm/branches/trunk-larour-NatGeoScience2016/src/c/modules/FloatingiceMeltingRatex/FloatingiceMeltingRatex.cpp

    r20020 r22481  
    2929                        MismipFloatingiceMeltingRatex(femmodel);
    3030                        break;
     31                case SpatialLinearFloatingMeltRateEnum:
     32                        if(VerboseSolution())_printf0_("        call Spatial Linear Floating melting rate module\n");
     33                        SpatialLinearFloatingiceMeltingRatex(femmodel);
     34                        break;
    3135                default:
    3236                        _error_("Basal forcing model "<<EnumToStringx(basalforcing_model)<<" not supported yet");
    … …  
    4347
    4448}/*}}}*/
     49void SpatialLinearFloatingiceMeltingRatex(FemModel* femmodel){/*{{{*/
     50
     51        for(int i=0;i<femmodel->elements->Size();i++){
     52                Element* element=xDynamicCast<Element*>(femmodel->elements->GetObjectByOffset(i));
     53                element->SpatialLinearFloatingiceMeltingRate();
     54        }
     55
     56}/*}}}*/
    4557void MismipFloatingiceMeltingRatex(FemModel* femmodel){/*{{{*/
    4658

  • issm/branches/trunk-larour-NatGeoScience2016/src/c/modules/FloatingiceMeltingRatex/FloatingiceMeltingRatex.h

    r19479 r22481  
    1111void FloatingiceMeltingRatex(FemModel* femmodel);
    1212void LinearFloatingiceMeltingRatex(FemModel* femmodel);
     13void SpatialLinearFloatingiceMeltingRatex(FemModel* femmodel);
    1314void MismipFloatingiceMeltingRatex(FemModel* femmodel);
    1415

  • issm/branches/trunk-larour-NatGeoScience2016/src/c/modules/ModelProcessorx/CreateParameters.cpp

    r22155 r22481  
    118118                        parameters->AddObject(iomodel->CopyConstantObject("md.basalforcings.upperwater_elevation",BasalforcingsUpperwaterElevationEnum));
    119119                        break;
     120                case SpatialLinearFloatingMeltRateEnum:
     121                        /*Nothing to add to parameters:*/
     122                        break;
    120123                case MismipFloatingMeltRateEnum:
    121124                        parameters->AddObject(iomodel->CopyConstantObject("md.basalforcings.meltrate_factor",BasalforcingsMeltrateFactorEnum));

  • issm/branches/trunk-larour-NatGeoScience2016/src/c/shared/Enum/EnumDefinitions.h

    r22248 r22481  
    7777        MismipFloatingMeltRateEnum,
    7878        MantlePlumeGeothermalFluxEnum,
     79        SpatialLinearFloatingMeltRateEnum,
    7980        BedEnum,
    8081        BaseEnum,

  • issm/branches/trunk-larour-NatGeoScience2016/src/c/shared/Enum/EnumToStringx.cpp

    r22248 r22481  
    8383                case MismipFloatingMeltRateEnum : return "MismipFloatingMeltRate";
    8484                case MantlePlumeGeothermalFluxEnum : return "MantlePlumeGeothermalFlux";
     85                case SpatialLinearFloatingMeltRateEnum : return "SpatialLinearFloatingMeltRate";
    8586                case BedEnum : return "Bed";
    8687                case BaseEnum : return "Base";

  • issm/branches/trunk-larour-NatGeoScience2016/src/c/shared/Enum/StringToEnumx.cpp

    r22248 r22481  
    8383              else if (strcmp(name,"MismipFloatingMeltRate")==0) return MismipFloatingMeltRateEnum;
    8484              else if (strcmp(name,"MantlePlumeGeothermalFlux")==0) return MantlePlumeGeothermalFluxEnum;
     85              else if (strcmp(name,"SpatialLinearFloatingMeltRate")==0) return SpatialLinearFloatingMeltRateEnum;
    8586              else if (strcmp(name,"Bed")==0) return BedEnum;
    8687              else if (strcmp(name,"Base")==0) return BaseEnum;
    … …  
    136137              else if (strcmp(name,"EffectivePressure")==0) return EffectivePressureEnum;
    137138              else if (strcmp(name,"EplHead")==0) return EplHeadEnum;
    138               else if (strcmp(name,"EplHeadOld")==0) return EplHeadOldEnum;
    139139         else stage=2;
    140140   }
    141141   if(stage==2){
    142               if (strcmp(name,"EplHeadSlopeX")==0) return EplHeadSlopeXEnum;
     142              if (strcmp(name,"EplHeadOld")==0) return EplHeadOldEnum;
     143              else if (strcmp(name,"EplHeadSlopeX")==0) return EplHeadSlopeXEnum;
    143144              else if (strcmp(name,"EplHeadSlopeY")==0) return EplHeadSlopeYEnum;
    144145              else if (strcmp(name,"EplZigZagCounter")==0) return EplZigZagCounterEnum;
    … …  
    259260              else if (strcmp(name,"CalvingMeltingrate")==0) return CalvingMeltingrateEnum;
    260261              else if (strcmp(name,"CalvingLevermann")==0) return CalvingLevermannEnum;
    261               else if (strcmp(name,"CalvingDev")==0) return CalvingDevEnum;
    262262         else stage=3;
    263263   }
    264264   if(stage==3){
    265               if (strcmp(name,"CalvingMinthickness")==0) return CalvingMinthicknessEnum;
     265              if (strcmp(name,"CalvingDev")==0) return CalvingDevEnum;
     266              else if (strcmp(name,"CalvingMinthickness")==0) return CalvingMinthicknessEnum;
    266267              else if (strcmp(name,"DefaultCalving")==0) return DefaultCalvingEnum;
    267268              else if (strcmp(name,"CalvinglevermannCoeff")==0) return CalvinglevermannCoeffEnum;
    … …  
    382383              else if (strcmp(name,"SmbSolution")==0) return SmbSolutionEnum;
    383384              else if (strcmp(name,"SmbNumRequestedOutputs")==0) return SmbNumRequestedOutputsEnum;
    384               else if (strcmp(name,"SmbRequestedOutputs")==0) return SmbRequestedOutputsEnum;
    385385         else stage=4;
    386386   }
    387387   if(stage==4){
    388               if (strcmp(name,"SmbIsInitialized")==0) return SmbIsInitializedEnum;
     388              if (strcmp(name,"SmbRequestedOutputs")==0) return SmbRequestedOutputsEnum;
     389              else if (strcmp(name,"SmbIsInitialized")==0) return SmbIsInitializedEnum;
    389390              else if (strcmp(name,"SmbDzini")==0) return SmbDziniEnum;
    390391              else if (strcmp(name,"SmbDini")==0) return SmbDiniEnum;
    … …  
    505506              else if (strcmp(name,"PressurePicard")==0) return PressurePicardEnum;
    506507              else if (strcmp(name,"AndroidFrictionCoefficient")==0) return AndroidFrictionCoefficientEnum;
    507               else if (strcmp(name,"ResetPenalties")==0) return ResetPenaltiesEnum;
    508508         else stage=5;
    509509   }
    510510   if(stage==5){
    511               if (strcmp(name,"SurfaceAbsVelMisfit")==0) return SurfaceAbsVelMisfitEnum;
     511              if (strcmp(name,"ResetPenalties")==0) return ResetPenaltiesEnum;
     512              else if (strcmp(name,"SurfaceAbsVelMisfit")==0) return SurfaceAbsVelMisfitEnum;
    512513              else if (strcmp(name,"SurfaceArea")==0) return SurfaceAreaEnum;
    513514              else if (strcmp(name,"SurfaceAverageVelMisfit")==0) return SurfaceAverageVelMisfitEnum;
    … …  
    628629              else if (strcmp(name,"Outputdefinition25")==0) return Outputdefinition25Enum;
    629630              else if (strcmp(name,"Outputdefinition26")==0) return Outputdefinition26Enum;
    630               else if (strcmp(name,"Outputdefinition27")==0) return Outputdefinition27Enum;
    631631         else stage=6;
    632632   }
    633633   if(stage==6){
    634               if (strcmp(name,"Outputdefinition28")==0) return Outputdefinition28Enum;
     634              if (strcmp(name,"Outputdefinition27")==0) return Outputdefinition27Enum;
     635              else if (strcmp(name,"Outputdefinition28")==0) return Outputdefinition28Enum;
    635636              else if (strcmp(name,"Outputdefinition29")==0) return Outputdefinition29Enum;
    636637              else if (strcmp(name,"Outputdefinition30")==0) return Outputdefinition30Enum;
    … …  
    751752              else if (strcmp(name,"ToolkitsOptionsAnalyses")==0) return ToolkitsOptionsAnalysesEnum;
    752753              else if (strcmp(name,"ToolkitsOptionsStrings")==0) return ToolkitsOptionsStringsEnum;
    753               else if (strcmp(name,"QmuErrName")==0) return QmuErrNameEnum;
    754754         else stage=7;
    755755   }
    756756   if(stage==7){
    757               if (strcmp(name,"QmuInName")==0) return QmuInNameEnum;
     757              if (strcmp(name,"QmuErrName")==0) return QmuErrNameEnum;
     758              else if (strcmp(name,"QmuInName")==0) return QmuInNameEnum;
    758759              else if (strcmp(name,"QmuOutName")==0) return QmuOutNameEnum;
    759760              else if (strcmp(name,"Regular")==0) return RegularEnum;
    … …  
    874875              else if (strcmp(name,"Matenhancedice")==0) return MatenhancediceEnum;
    875876              else if (strcmp(name,"Matestar")==0) return MatestarEnum;
    876               else if (strcmp(name,"Matpar")==0) return MatparEnum;
    877877         else stage=8;
    878878   }
    879879   if(stage==8){
    880               if (strcmp(name,"Node")==0) return NodeEnum;
     880              if (strcmp(name,"Matpar")==0) return MatparEnum;
     881              else if (strcmp(name,"Node")==0) return NodeEnum;
    881882              else if (strcmp(name,"Numericalflux")==0) return NumericalfluxEnum;
    882883              else if (strcmp(name,"Neumannflux")==0) return NeumannfluxEnum;
    … …  
    997998              else if (strcmp(name,"MaxAbsVz")==0) return MaxAbsVzEnum;
    998999              else if (strcmp(name,"FloatingArea")==0) return FloatingAreaEnum;
    999               else if (strcmp(name,"GroundedArea")==0) return GroundedAreaEnum;
    10001000         else stage=9;
    10011001   }
    10021002   if(stage==9){
    1003               if (strcmp(name,"IceMass")==0) return IceMassEnum;
     1003              if (strcmp(name,"GroundedArea")==0) return GroundedAreaEnum;
     1004              else if (strcmp(name,"IceMass")==0) return IceMassEnum;
    10041005              else if (strcmp(name,"IceVolume")==0) return IceVolumeEnum;
    10051006              else if (strcmp(name,"IceVolumeAboveFloatation")==0) return IceVolumeAboveFloatationEnum;

  • issm/branches/trunk-larour-NatGeoScience2016/src/c/shared/io/Marshalling/IoCodeConversions.cpp

    r21759 r22481  
    9292                case 3: return MismipFloatingMeltRateEnum;
    9393                case 4: return MantlePlumeGeothermalFluxEnum;
     94                case 5: return SpatialLinearFloatingMeltRateEnum;
    9495                default: _error_("Marshalled Basal Forcings code \""<<enum_in<<"\" not supported yet");
    9596        }

  • issm/branches/trunk-larour-NatGeoScience2016/src/m/classes/model.m

    r22091 r22481  
    328328
    329329                end % }}}
    330                 function md2 = extract(md,area) % {{{
     330                function md2 = extract(md,area,varargin) % {{{
    331331                        %extract - extract a model according to an Argus contour or flag list
    332332                        %
    … …  
    350350                        md1=md;
    351351
     352                        %recover optoins:
     353                        options=pairoptions(varargin{:});
     354
    352355                        %some checks
    353                         if ((nargin~=2) | (nargout~=1)),
     356                        if ((nargin<2) | (nargout~=1)),
    354357                                help extract
    355358                                error('extract error message: bad usage');
    … …  
    363366
    364367                        %kick out all elements with 3 dirichlets
    365                         spc_elem=find(~flag_elem);
    366                         spc_node=sort(unique(md1.mesh.elements(spc_elem,:)));
    367                         flag=ones(md1.mesh.numberofvertices,1);
    368                         flag(spc_node)=0;
    369                         pos=find(sum(flag(md1.mesh.elements),2)==0);
    370                         flag_elem(pos)=0;
     368                        if getfieldvalue(options,'spccheck',1)
     369                                spc_elem=find(~flag_elem);
     370                                spc_node=sort(unique(md1.mesh.elements(spc_elem,:)));
     371                                flag=ones(md1.mesh.numberofvertices,1);
     372                                flag(spc_node)=0;
     373                                pos=find(sum(flag(md1.mesh.elements),2)==0);
     374                                flag_elem(pos)=0;
     375                        end
    371376
    372377                        %extracted elements and nodes lists
    … …  
    13891394
    13901395                end % }}}
    1391                 function savemodeljs(md,modelname,websiteroot) % {{{
    1392 
     1396                function savemodeljs(md,modelname,websiteroot,varargin) % {{{
     1397               
    13931398                        %the goal of this routine is to save the model as a javascript array that can be included in any html
    13941399                        %file:
    13951400
     1401                        options=pairoptions(varargin{:});
     1402                        optimization=getfieldvalue(options,'optimize',0);
     1403
     1404                       
    13961405                        %disp:
    13971406                        disp(['saving model ''' modelname ''' in file ' websiteroot '/js/' modelname '.js']);
    … …  
    14121421                                end
    14131422
     1423                                %some optimization:
     1424                                if optimization==1,
     1425                                        %optimize for plotting only:
     1426                                        if ~ismember(field,{'geometry','mesh','mask'}),
     1427                                                continue;
     1428                                        end
     1429                                end
     1430
    14141431                                %Check that current field is an object
    14151432                                if ~isobject(md.(field))

  • issm/branches/trunk-larour-NatGeoScience2016/src/m/classes/sealevelmodel.m

    r22244 r22481  
    133133                        end
    134134                end % }}}
     135                function listcaps(self) % {{{
     136                        for  i=1:length(self.icecaps),
     137                                disp(sprintf('%i: %s',i,self.icecaps{i}.miscellaneous.name));
     138                        end
     139                end % }}}
     140                function viscousiterations(self) % {{{
     141                        for  i=1:length(self.icecaps),
     142                                ic=self.icecaps{i};
     143                                mvi=ic.results.TransientSolution(1).StressbalanceConvergenceNumSteps;
     144                                for j=2:length(ic.results.TransientSolution)-1,
     145                                        mvi=max(mvi,ic.results.TransientSolution(j).StressbalanceConvergenceNumSteps);
     146                                end
     147                                disp(sprintf('%i, %s: %i',i,self.icecaps{i}.miscellaneous.name,mvi));
     148                        end
     149                end % }}}
     150                function maxtimestep(self) % {{{
     151                        for  i=1:length(self.icecaps),
     152                                ic=self.icecaps{i};
     153                                mvi=length(ic.results.TransientSolution);
     154                                timei=ic.results.TransientSolution(end).time;
     155                                disp(sprintf('%i, %s: %i/%g',i,self.icecaps{i}.miscellaneous.name,mvi,timei));
     156                        end
     157                        mvi=length(self.earth.results.TransientSolution);
     158                        timei=self.earth.results.TransientSolution(end).time;
     159                        disp(sprintf('Earth: %i/%g',mvi,timei));
     160                end % }}}
     161                function self=homogeneize(self) % {{{
     162                        mintimestep=Inf;
     163                        for  i=1:length(self.icecaps),
     164                                ic=self.icecaps{i};
     165                                mintimestep=min(mintimestep, length(ic.results.TransientSolution));
     166                        end
     167                        mintimestep=min(mintimestep, length(self.earth.results.TransientSolution));
     168                       
     169                        for  i=1:length(self.icecaps),
     170                                ic=self.icecaps{i};
     171                                ic.results.TransientSolution=ic.results.TransientSolution(1:mintimestep);
     172                                self.icecaps{i}=ic;
     173                        end
     174                        ic=self.earth;
     175                        ic.results.TransientSolution=ic.results.TransientSolution(1:mintimestep);
     176                        self.earth=ic;
     177                end % }}}
    135178        end
    136179end

  • issm/branches/trunk-larour-NatGeoScience2016/src/m/exp/contourlevelzero.m

    r17686 r22481  
    1616x=md.mesh.x;
    1717y=md.mesh.y;
     18if isfield(md.mesh,'z'),
     19        z=md.mesh.z;
     20else
     21        z=zeros(md.mesh.numberofvertices,1);
     22end
    1823index=md.mesh.elements;
    1924
    … …  
    8287y1=zeros(numelems,1);
    8388y2=zeros(numelems,1);
     89z1=zeros(numelems,1);
     90z2=zeros(numelems,1);
     91
    8492edge_l=zeros(numelems,2);
    8593
    … …  
    96104                y1(j)=y(Seg1(poselem(j),1))+weight1*(y(Seg1(poselem(j),2))-y(Seg1(poselem(j),1)));
    97105                y2(j)=y(Seg2(poselem(j),1))+weight2*(y(Seg2(poselem(j),2))-y(Seg2(poselem(j),1)));
     106                z1(j)=z(Seg1(poselem(j),1))+weight1*(z(Seg1(poselem(j),2))-z(Seg1(poselem(j),1)));
     107                z2(j)=z(Seg2(poselem(j),1))+weight2*(z(Seg2(poselem(j),2))-z(Seg2(poselem(j),1)));
     108
    98109                edge_l(j,1)=Seg1_num(poselem(j));
    99110                edge_l(j,2)=Seg2_num(poselem(j));
    … …  
    105116                y1(j)=y(Seg1(poselem(j),1))+weight1*(y(Seg1(poselem(j),2))-y(Seg1(poselem(j),1)));
    106117                y2(j)=y(Seg3(poselem(j),1))+weight3*(y(Seg3(poselem(j),2))-y(Seg3(poselem(j),1)));
     118                z1(j)=z(Seg1(poselem(j),1))+weight1*(z(Seg1(poselem(j),2))-z(Seg1(poselem(j),1)));
     119                z2(j)=z(Seg3(poselem(j),1))+weight3*(z(Seg3(poselem(j),2))-z(Seg3(poselem(j),1)));
     120
    107121                edge_l(j,1)=Seg1_num(poselem(j));
    108122                edge_l(j,2)=Seg3_num(poselem(j));
    … …  
    114128                y1(j)=y(Seg2(poselem(j),1))+weight2*(y(Seg2(poselem(j),2))-y(Seg2(poselem(j),1)));
    115129                y2(j)=y(Seg3(poselem(j),1))+weight3*(y(Seg3(poselem(j),2))-y(Seg3(poselem(j),1)));
     130                z1(j)=z(Seg2(poselem(j),1))+weight2*(z(Seg2(poselem(j),2))-z(Seg2(poselem(j),1)));
     131                z2(j)=z(Seg3(poselem(j),1))+weight3*(z(Seg3(poselem(j),2))-z(Seg3(poselem(j),1)));
     132
    116133                edge_l(j,1)=Seg2_num(poselem(j));
    117134                edge_l(j,2)=Seg3_num(poselem(j));
    … …  
    130147        %take the right edge of the second segment and connect it to the next segments if any
    131148        e1=edge_l(1,1);   e2=edge_l(1,2);
    132         xc=[x1(1);x2(1)]; yc=[y1(1);y2(1)];
     149        xc=[x1(1);x2(1)]; yc=[y1(1);y2(1)]; zc=[z1(1);z2(1)];
     150
    133151
    134152        %erase the lines corresponding to this edge
    … …  
    136154        x1(1)=[]; x2(1)=[];
    137155        y1(1)=[]; y2(1)=[];
     156        z1(1)=[]; z2(1)=[];
    138157
    139158        [ro1,co1]=find(edge_l==e1);
    … …  
    142161
    143162                if co1==1,
    144                         xc=[x2(ro1);xc]; yc=[y2(ro1);yc];
     163                        xc=[x2(ro1);xc]; yc=[y2(ro1);yc];zc=[z2(ro1);zc];
    145164
    146165                        %next edge:
    … …  
    148167
    149168                else
    150                         xc=[x1(ro1);xc]; yc=[y1(ro1);yc];
     169                        xc=[x1(ro1);xc]; yc=[y1(ro1);yc];zc=[z1(ro1);zc];
    151170
    152171                        %next edge:
    … …  
    158177                x1(ro1)=[]; x2(ro1)=[];
    159178                y1(ro1)=[]; y2(ro1)=[];
     179                z1(ro1)=[]; z2(ro1)=[];
    160180
    161181                %next connection
    … …  
    169189
    170190                if co2==1,
    171                         xc=[xc;x2(ro2)]; yc=[yc;y2(ro2)];
     191                        xc=[xc;x2(ro2)]; yc=[yc;y2(ro2)];zc=[zc;z2(ro2)];
    172192
    173193                        %next edge:
    174194                        e2=edge_l(ro2,2);
    175195                else
    176                         xc=[xc;x1(ro2)]; yc=[yc;y1(ro2)];
     196                        xc=[xc;x1(ro2)]; yc=[yc;y1(ro2)]; zc=[zc;z1(ro2)];
    177197
    178198                        %next edge:
    … …  
    184204                x1(ro2)=[]; x2(ro2)=[];
    185205                y1(ro2)=[]; y2(ro2)=[];
     206                z1(ro2)=[]; z2(ro2)=[];
    186207
    187208                %next connection
    … …  
    192213        contours(end+1).x=xc;
    193214        contours(end).y=yc;
     215        contours(end).z=zc;
    194216        contours(end).name='';
    195217        contours(end).nods=length(xc);

  • issm/branches/trunk-larour-NatGeoScience2016/src/m/mesh/meshintersect.m

    r22217 r22481  
    2828                        end
    2929                        error('one or more vertices on the global mesh were duplicated!');
     30                elseif isempty(s),
     31                        plot(distance)
     32                        min(distance)
     33                        error('could not find vertices in common, relax tolerance?');
    3034                else
    3135                        indices(i)=s;

  • issm/branches/trunk-larour-NatGeoScience2016/src/m/mesh/modelmerge2d.m

    r22217 r22481  
    1010        tolerance=getfieldvalue(options,'tolerance',10^-5);
    1111       
    12         md=model();
     12        md=md1; %by default, we transfer all the settings from md1 to md.
    1313
    1414        %first ,copy md1 mesh into md.mesh to initialize:
    … …  
    8383        md.mesh.vertexonboundary(md.mesh.segments(:,1:2))=1;
    8484
     85        if getfieldvalue(options,'full',0),
     86                %we are asked to merge the classes fields too. We need have vertex and element mappings first:
     87                %vertex intersections:
     88                md.mesh.extractedvertices={meshintersect(x,y,md1.mesh.x,md1.mesh.y,'tolerance',1e-5), meshintersect(x,y,md2.mesh.x,md2.mesh.y,'tolerance',1e-5)};
     89                %element intersections:
     90                xe=x(md.mesh.elements)*[1;1;1]/3; ye=y(md.mesh.elements)*[1;1;1]/3;
     91                x1e=md1.mesh.x(md1.mesh.elements)*[1;1;1]/3; y1e=md1.mesh.y(md1.mesh.elements)*[1;1;1]/3;
     92                x2e=md2.mesh.x(md2.mesh.elements)*[1;1;1]/3; y2e=md2.mesh.y(md2.mesh.elements)*[1;1;1]/3;
     93                md.mesh.extractedelements= {meshintersect(xe,ye,x1e,y1e,'tolerance',1e-5) , meshintersect(xe,ye,x2e,y2e,'tolerance',1e-5)};
     94
     95                %now we can go through classes and transfer.
     96                md=transfer_fields(md,md1,md2,'geometry',{'thickness','surface','bed','base'});
     97                md=transfer_fields(md,md1,md2,'mask',{'groundedice_levelset','ice_levelset','ocean_levelset','land_levelset','glacier_levelset'});
     98                md=transfer_fields(md,md1,md2,'smb',{'mass_balance'});
     99                if strcmpi(class(md1.basalforcings),'linearbasalforcings'),
     100                        md=transfer_fields(md,md1,md2,'basalforcings',{'groundedice_melting_rate','geothermalflux'});
     101                else
     102                        md=transfer_fields(md,md1,md2,'basalforcings',{'groundedice_melting_rate','deepwater_melting_rate','deepwater_elevation','upperwater_elevation','geothermalflux'});
     103                end
     104                md=transfer_fields(md,md1,md2,'materials',{'rheology_B','rheology_n'});
     105                md=transfer_fields(md,md1,md2,'friction',{'coefficient','p','q'});
     106                md=transfer_fields(md,md1,md2,'flowequation',{'vertex_equation','element_equation','borderSSA','borderFS','borderHO'});
     107                md=transfer_fields(md,md1,md2,'initialization',{'vx','vy','vz','vel','pressure','temperature'});
     108                md=transfer_fields(md,md1,md2,'slr',{'deltathickness','sealevel','spcthickness','steric_rate'});
     109                md=transfer_fields(md,md1,md2,'masstransport',{'spcthickness'});
     110                md=transfer_fields(md,md1,md2,'thermal',{'spctemperature'});
     111                md=transfer_fields(md,md1,md2,'inversion',{'min_parameters','max_parameters','vx_obs','vy_obs','vz_obs'});
     112                md.inversion.cost_functions_coefficients=zeros(md.mesh.numberofvertices,3);
     113                md.inversion.cost_functions_coefficients(md.mesh.extractedvertices{1},:)=md1.inversion.cost_functions_coefficients;
     114                md.inversion.cost_functions_coefficients(md.mesh.extractedvertices{2},:)=md2.inversion.cost_functions_coefficients;
     115
     116                %boundary conditions:
     117                md=transfer_fields(md,md1,md2,'stressbalance',{'spcvx','spcvy','spcvz'});
     118                md.stressbalance.loadingforce=zeros(md.mesh.numberofvertices,3);
     119                md.stressbalance.referential=NaN*ones(md.mesh.numberofvertices,6);
     120                bound1=zeros(md.mesh.numberofvertices,1); bound1(md.mesh.extractedvertices{1})=md1.mesh.vertexonboundary;
     121                bound2=zeros(md.mesh.numberofvertices,1); bound2(md.mesh.extractedvertices{2})=md2.mesh.vertexonboundary;
     122                boundary=bound1 & bound2;
     123                pos=find(boundary); md.stressbalance.spcvx(pos)=NaN; md.stressbalance.spcvy(pos)=NaN; md.stressbalance.spcvz(pos)=NaN;
     124        end
     125
     126       
    85127        %some checks:
    86128        if max(md.mesh.elements)>md.mesh.numberofvertices,
    87129                error('issue in modelmerge, one of the element ids is > number of vertices!');
    88130        end
     131
     132end %end of function
     133
     134function prop=transfer_vertices(md,md1,md2,field1,field2) % {{{
     135        f1=getfield(md1,field1); f2=getfield(f1,field2);
     136        if length(f2)==md1.mesh.numberofvertices,
     137                prop=zeros(md.mesh.numberofvertices,1);
     138                prop(md.mesh.extractedvertices{1})=f2;
     139                f1=getfield(md2,field1); f2=getfield(f1,field2); prop(md.mesh.extractedvertices{2})=f2;
     140        else
     141                prop=zeros(md.mesh.numberofvertices+1,1);  prop(end)=f2(end);
     142                prop(md.mesh.extractedvertices{1})=f2(1:end-1);
     143                f1=getfield(md2,field1); f2=getfield(f1,field2); prop(md.mesh.extractedvertices{2})=f2(1:end-1);
     144                prop=zeros(md.mesh.numberofvertices+1,1);
     145        end
     146       
     147       
     148end %end of function %}}}
     149function prop=transfer_elements(md,md1,md2,field1,field2) % {{{
     150        prop=zeros(md.mesh.numberofelements,1);
     151        f1=getfield(md1,field1); f2=getfield(f1,field2); prop(md.mesh.extractedelements{1})=f2;
     152        f1=getfield(md2,field1); f2=getfield(f1,field2); prop(md.mesh.extractedelements{2})=f2;
     153
     154end %end of function %}}}
     155function md=transfer_fields(md,md1,md2,classname,classfields) % {{{
     156
     157        for i=1:length(classfields),
     158                field1=eval(['md1.' classname '.' classfields{i}]);
     159                if length(field1)==md1.mesh.numberofvertices | length(field1)==md1.mesh.numberofvertices+1,
     160                        eval(['md.' classname '.' classfields{i} '=transfer_vertices(md,md1,md2,''' classname ''',''' classfields{i} ''');']);
     161                else
     162                        eval(['md.' classname '.' classfields{i} '=transfer_elements(md,md1,md2,''' classname ''',''' classfields{i} ''');']);
     163                end
     164        end
     165
     166end %end of function %}}}

  • issm/branches/trunk-larour-NatGeoScience2016/src/m/plot/applyoptions.m

    r22218 r22481  
    265265        end
    266266end
     267
     268%shpdisp
     269if exist(options,'shpdisp'),
     270        filename=(getfieldvalue(options,'shpdisp'));
     271        style=(getfieldvalue(options,'shpstyle',{'r.-'}));
     272        linewidth=(getfieldvalue(options,'linewidth',1));
     273        for i=1:length(getfieldvalue(options,'shpdisp')),
     274                filenamei=filename{i};
     275                stylei=style{i};
     276                if length(linewidth)==1,
     277                        linewidthi=linewidth;
     278                else
     279                        linewidthi=linewidth{i};
     280                end
     281                %shpdisp(filenamei,'linestyle',stylei,'linewidth',linewidthi,'multiplier',getfieldvalue(options,'unit',1));
     282                shpdisp(filenamei,1,stylei,linewidthi,getfieldvalue(options,'unit',1));
     283        end
     284end
     285
     286
    267287
    268288%text (default value is empty, not NaN...)

  • issm/branches/trunk-larour-NatGeoScience2016/src/m/plot/googlemaps.m

    r21759 r22481  
    121121
    122122%Write image
    123 imwrite(final,'temp.png','png')
     123imwrite(final,'temp.png','png');
    124124[ulmx ulmy]=ll2mercator(ullat,ullon);
    125125[lrmx lrmy]=ll2mercator(lrlat,lrlon);

  • issm/branches/trunk-larour-NatGeoScience2016/src/m/plot/plot_contour.m

    r22218 r22481  
    2121if datatype==1,
    2222        %elements -> take average
    23         data=averaging(md,data,0)
     23        data=averaging(md,data,0);
    2424elseif datatype==2,
    2525        %nodes -> do nothing

  • issm/branches/trunk-larour-NatGeoScience2016/src/m/plot/plot_riftrelvel.m

    r13870 r22481  
    110110hold off
    111111
     112%xlim and ylim:
     113if ~exist(options,'xlim') options=addfielddefault(options,'xlim',[min(md.mesh.x) max(md.mesh.x)]); end
     114if ~exist(options,'ylim') options=addfielddefault(options,'ylim',[min(md.mesh.y) max(md.mesh.y)]); end
     115
    112116%apply options
    113117quiver_colorbar(quivers,options);

  • issm/branches/trunk-larour-NatGeoScience2016/src/m/plot/plot_rifts.m

    r19064 r22481  
    2424                        normal(1)=penaltypairs(j,5);
    2525                        normal(2)=penaltypairs(j,6);
    26                         x(penaltypairs(j,1))=x(penaltypairs(j,1))-normal(1)*offset;
    27                         y(penaltypairs(j,1))=y(penaltypairs(j,1))-normal(2)*offset;
     26                        x(penaltypairs(j,1))=x(penaltypairs(j,1))-normal(1)*offset/2;
     27                        y(penaltypairs(j,1))=y(penaltypairs(j,1))-normal(2)*offset/2;
     28                        x(penaltypairs(j,2))=x(penaltypairs(j,2))+normal(1)*offset/2;
     29                        y(penaltypairs(j,2))=y(penaltypairs(j,2))+normal(2)*offset/2;
     30
    2831                end
    2932                if length(md.rifts.riftstruct(i).tips)==3,

  • issm/branches/trunk-larour-NatGeoScience2016/src/m/plot/plot_riftvel.m

    r13730 r22481  
    101101hold off
    102102
     103%xlim and ylim:
     104if ~exist(options,'xlim') options=addfielddefault(options,'xlim',[min(md.mesh.x) max(md.mesh.x)]); end
     105if ~exist(options,'ylim') options=addfielddefault(options,'ylim',[min(md.mesh.y) max(md.mesh.y)]); end
     106
    103107%apply options
    104108quiver_colorbar(quivers,options);

  • issm/branches/trunk-larour-NatGeoScience2016/src/m/plot/showregion.m

    r18558 r22481  
    2929mdy=[min(md.mesh.y) max(md.mesh.y)];
    3030
    31 line(A.x,A.y,'color','b');
    3231patch([Ax(1)  Ax(2)  Ax(2)  Ax(1) Ax(1)],[Ay(1)  Ay(1)  Ay(2)  Ay(2) Ay(1)],[1 1 1],'EdgeColor',[0 0 0],'LineWidth',1,'FaceLighting','none')
    3332patch( [mdx(1) mdx(2) mdx(2) mdx(1)],[mdy(1) mdy(1) mdy(2) mdy(2)],ones(4,1),'EdgeColor',[0 0 0],'FaceColor','r','FaceAlpha',0.5)
     33line(A.x,A.y,'color','k');
    3434colorbar('off');

  • issm/branches/trunk-larour-NatGeoScience2016/src/m/solve/loadresultsfromcluster.m

    r19335 r22481  
    1 function md=loadresultsfromcluster(md,runtimename)
     1function md=loadresultsfromcluster(md,varargin)
    22%LOADRESULTSFROMCLUSTER - load results of solution sequence from cluster
    33%
    44%   Usage:
    5 %      md=loadresultsfromcluster(md,runtimename);
     5%      md=loadresultsfromcluster(md,varargin);
     6%      options include runtimename and nolog.
     7
     8%process options:
     9options=pairoptions(varargin{:});
     10nolog=getfieldvalue(options,'nolog',0);
    611
    712%retrieve cluster, to be able to call its methods
    813cluster=md.cluster;
    914
    10 if nargin==2,
    11         md.private.runtimename=runtimename;
     15if exist(options,'runtimename'),
     16        md.private.runtimename=getfieldvalue(optoins,'runtimename');
    1217end
    1318
    1419%Download outputs from the cluster
    15 filelist={[md.miscellaneous.name '.outlog'],[md.miscellaneous.name '.errlog']};
     20if ~nolog,
     21        filelist={[md.miscellaneous.name '.outlog'],[md.miscellaneous.name '.errlog']};
     22else
     23        filelist={};
     24end
    1625if md.qmu.isdakota,
    1726        filelist{end+1}=[md.miscellaneous.name '.qmu.err'];
    … …  
    3544        delete([['qmu' num2str(feature('GetPid')) '/']  md.miscellaneous.name '.errlog']);
    3645else
    37         delete([md.miscellaneous.name '.outlog']);
    38         delete([md.miscellaneous.name '.errlog']);
     46        if ~nolog,
     47                delete([md.miscellaneous.name '.outlog']);
     48                delete([md.miscellaneous.name '.errlog']);
     49        end
    3950        delete([md.miscellaneous.name '.outbin']);
    4051        if exist([md.private.runtimename '.tar.gz']) & ~ispc(),

  • issm/branches/trunk-larour-NatGeoScience2016/src/m/solve/loadresultsfromclusterslm.m

    r21323 r22481  
    1010        %go through icecaps and download results
    1111        for i=1:length(slm.icecaps),
    12                 slm.icecaps{i}=loadresultsfromcluster(slm.icecaps{i});
     12                slm.icecaps{i}=loadresultsfromcluster(slm.icecaps{i},'nolog',1);
    1313        end
    14         slm.earth=loadresultsfromcluster(slm.earth);
     14        slm.earth=loadresultsfromcluster(slm.earth,'nolog',1);
Note: See TracChangeset for help on using the changeset viewer.