Index: /issm/trunk-jpl/src/c/classes/Elements/Element.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Element.cpp (revision 21807) +++ /issm/trunk-jpl/src/c/classes/Elements/Element.cpp (revision 21808) @@ -1341,4 +1341,17 @@ } /*}}}*/ +IssmDouble Element::GroundedArea(IssmDouble* mask){/*{{{*/ + + /*Retrieve values of the mask defining the element: */ + for(int i=0;iGetNumberOfVertices();i++){ + if(mask[this->vertices[i]->Sid()]<=0.){ + return 0.; + } + } + + /*Return: */ + return this->GroundedArea(); +} +/*}}}*/ bool Element::HasNodeOnBase(){/*{{{*/ return (this->inputs->Max(MeshVertexonbaseEnum)>0.); @@ -1347,4 +1360,55 @@ return (this->inputs->Max(MeshVertexonsurfaceEnum)>0.); }/*}}}*/ +IssmDouble Element::IceMass(void){/*{{{*/ + + IssmDouble rho_ice; + + if(!IsIceInElement())return 0.; //do not contribute to the volume of the ice! + + /*recover ice density: */ + rho_ice=matpar->GetMaterialParameter(MaterialsRhoIceEnum); + + return rho_ice*this->IceVolume(); +} +/*}}}*/ +IssmDouble Element::IceMass(IssmDouble* mask){/*{{{*/ + + /*Retrieve values of the mask defining the element: */ + for(int i=0;iGetNumberOfVertices();i++){ + if(mask[this->vertices[i]->Sid()]<=0.){ + return 0.; + } + } + + /*Return: */ + return this->IceMass(); +} +/*}}}*/ +IssmDouble Element::IceVolume(IssmDouble* mask){/*{{{*/ + + /*Retrieve values of the mask defining the element: */ + for(int i=0;iGetNumberOfVertices();i++){ + if(mask[this->vertices[i]->Sid()]<=0.){ + return 0.; + } + } + + /*Return: */ + return this->IceVolume(); +} +/*}}}*/ +IssmDouble Element::IceVolumeAboveFloatation(IssmDouble* mask){/*{{{*/ + + /*Retrieve values of the mask defining the element: */ + for(int i=0;iGetNumberOfVertices();i++){ + if(mask[this->vertices[i]->Sid()]<=0.){ + return 0.; + } + } + + /*Return: */ + return this->IceVolumeAboveFloatation(); +} +/*}}}*/ int Element::Id(){/*{{{*/ @@ -2952,4 +3016,43 @@ return this->matpar->TMeltingPoint(pressure); }/*}}}*/ +IssmDouble Element::TotalFloatingBmb(IssmDouble* mask){/*{{{*/ + + /*Retrieve values of the mask defining the element: */ + for(int i=0;iGetNumberOfVertices();i++){ + if(mask[this->vertices[i]->Sid()]<=0.){ + return 0.; + } + } + + /*Return: */ + return this->TotalFloatingBmb(); +} +/*}}}*/ +IssmDouble Element::TotalGroundedBmb(IssmDouble* mask){/*{{{*/ + + /*Retrieve values of the mask defining the element: */ + for(int i=0;iGetNumberOfVertices();i++){ + if(mask[this->vertices[i]->Sid()]<=0.){ + return 0.; + } + } + + /*Return: */ + return this->TotalGroundedBmb(); +} +/*}}}*/ +IssmDouble Element::TotalSmb(IssmDouble* mask){/*{{{*/ + + /*Retrieve values of the mask defining the element: */ + for(int i=0;iGetNumberOfVertices();i++){ + if(mask[this->vertices[i]->Sid()]<=0.){ + return 0.; + } + } + + /*Return: */ + return this->TotalSmb(); +} +/*}}}*/ void Element::TransformInvStiffnessMatrixCoord(ElementMatrix* Ke,int transformenum){/*{{{*/ Index: /issm/trunk-jpl/src/c/classes/Elements/Element.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Element.h (revision 21807) +++ /issm/trunk-jpl/src/c/classes/Elements/Element.h (revision 21808) @@ -109,6 +109,11 @@ IssmDouble GetZcoord(IssmDouble* xyz_list,Gauss* gauss); void GradientIndexing(int* indexing,int control_index,bool onsid=false); + IssmDouble GroundedArea(IssmDouble* mask); bool HasNodeOnBase(); bool HasNodeOnSurface(); + IssmDouble IceMass(); + IssmDouble IceMass(IssmDouble* mask); + IssmDouble IceVolume(IssmDouble* mask); + IssmDouble IceVolumeAboveFloatation(IssmDouble* mask); int Id(); void InputChangeName(int enum_type,int enum_type_old); @@ -150,4 +155,7 @@ void ThermalToEnthalpy(IssmDouble* penthalpy,IssmDouble temperature,IssmDouble waterfraction,IssmDouble pressure); IssmDouble TMeltingPoint(IssmDouble pressure); + IssmDouble TotalFloatingBmb(IssmDouble* mask); + IssmDouble TotalGroundedBmb(IssmDouble* mask); + IssmDouble TotalSmb(IssmDouble* mask); void TransformInvStiffnessMatrixCoord(ElementMatrix* Ke,int cs_enum); void TransformInvStiffnessMatrixCoord(ElementMatrix* Ke,Node** nodes,int numnodes,int cs_enum); @@ -212,5 +220,4 @@ virtual void GetVerticesCoordinatesTop(IssmDouble** xyz_list)=0; virtual IssmDouble GroundedArea(void)=0; - virtual IssmDouble IceMass(void)=0; virtual IssmDouble IceVolume(void)=0; virtual IssmDouble IceVolumeAboveFloatation(void)=0; Index: /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp (revision 21807) +++ /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp (revision 21808) @@ -1223,16 +1223,4 @@ /*Clean up and return*/ return phi*base_area; -} -/*}}}*/ -IssmDouble Penta::IceMass(void){/*{{{*/ - - IssmDouble rho_ice; - - if(!IsIceInElement())return 0.; //do not contribute to the volume of the ice! - - /*recover ice density: */ - rho_ice=matpar->GetMaterialParameter(MaterialsRhoIceEnum); - - return rho_ice*this->IceVolume(); } /*}}}*/ Index: /issm/trunk-jpl/src/c/classes/Elements/Penta.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Penta.h (revision 21807) +++ /issm/trunk-jpl/src/c/classes/Elements/Penta.h (revision 21808) @@ -92,5 +92,4 @@ void GetVerticesCoordinatesTop(IssmDouble** pxyz_list); IssmDouble GroundedArea(void); - IssmDouble IceMass(void); IssmDouble IceVolume(void); IssmDouble IceVolumeAboveFloatation(void); Index: /issm/trunk-jpl/src/c/classes/Elements/Seg.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Seg.h (revision 21807) +++ /issm/trunk-jpl/src/c/classes/Elements/Seg.h (revision 21808) @@ -79,5 +79,4 @@ void GetVerticesCoordinatesTop(IssmDouble** pxyz_list){_error_("not implemented yet");}; IssmDouble GroundedArea(void){_error_("not implemented yet");}; - IssmDouble IceMass(void){_error_("not implemented yet");}; IssmDouble IceVolume(void){_error_("not implemented yet");}; IssmDouble IceVolumeAboveFloatation(void){_error_("not implemented yet");}; Index: /issm/trunk-jpl/src/c/classes/Elements/Tetra.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Tetra.h (revision 21807) +++ /issm/trunk-jpl/src/c/classes/Elements/Tetra.h (revision 21808) @@ -87,5 +87,4 @@ bool HasFaceOnBase(); bool HasFaceOnSurface(); - IssmDouble IceMass(void){_error_("not implemented yet");}; IssmDouble IceVolume(void){_error_("not implemented yet");}; IssmDouble IceVolumeAboveFloatation(void){_error_("not implemented yet");}; Index: /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp (revision 21807) +++ /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp (revision 21808) @@ -1588,16 +1588,4 @@ return false; } -} -/*}}}*/ -IssmDouble Tria::IceMass(void){/*{{{*/ - - IssmDouble rho_ice; - - if(!IsIceInElement())return 0.; //do not contribute to the volume of the ice! - - /*recover ice density: */ - rho_ice=matpar->GetMaterialParameter(MaterialsRhoIceEnum); - - return rho_ice*this->IceVolume(); } /*}}}*/ Index: /issm/trunk-jpl/src/c/classes/Elements/Tria.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Tria.h (revision 21807) +++ /issm/trunk-jpl/src/c/classes/Elements/Tria.h (revision 21808) @@ -91,5 +91,4 @@ bool HasEdgeOnBase(); bool HasEdgeOnSurface(); - IssmDouble IceMass(void); IssmDouble IceVolume(void); IssmDouble IceVolumeAboveFloatation(void); Index: /issm/trunk-jpl/src/c/classes/Massfluxatgate.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Massfluxatgate.h (revision 21807) +++ /issm/trunk-jpl/src/c/classes/Massfluxatgate.h (revision 21808) @@ -84,9 +84,11 @@ /*}}}*/ ~Massfluxatgate(){/*{{{*/ - xDelete(this->x1); - xDelete(this->y1); - xDelete(this->x2); - xDelete(this->y2); - xDelete(this->elements); + if(this->numsegments){ + xDelete(this->x1); + xDelete(this->y1); + xDelete(this->x2); + xDelete(this->y2); + xDelete(this->elements); + } xDelete(this->name); } Index: /issm/trunk-jpl/src/c/classes/Regionaloutput.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Regionaloutput.h (revision 21808) +++ /issm/trunk-jpl/src/c/classes/Regionaloutput.h (revision 21808) @@ -0,0 +1,148 @@ +/*!\file Regionaloutput.h + * \brief: header file for Regionaloutput object + */ + +#ifndef _REGIONALOUTPUT_H_ +#define _REGIONALOUTPUT_H_ + +/*Headers:*/ +/*{{{*/ +#include "./Definition.h" +#include "../datastructures/datastructures.h" +#include "./Elements/Element.h" +#include "./Elements/Elements.h" +#include "./FemModel.h" +#include "../classes/Params/Parameters.h" + +/*}}}*/ + +class Regionaloutput: public Object, public Definition{ + + public: + + int definitionenum; + char* outputname; + char* name; + IssmDouble* mask; + int M; + + /*Regionalicevolume: constructors, destructors :*/ + Regionaloutput(){/*{{{*/ + + this->definitionenum = -1; + this->outputname = NULL; + this->name = NULL; + this->mask=NULL; + this->M=0; + + } + /*}}}*/ + Regionaloutput(char* in_name, int in_definitionenum, char* in_outputname, IssmDouble* maskin, int Min){ /*{{{*/ + + this->definitionenum=in_definitionenum; + this->outputname = xNew(strlen(in_outputname)+1); + xMemCpy(this->outputname,in_outputname,strlen(in_outputname)+1); + this->name = xNew(strlen(in_name)+1); + xMemCpy(this->name,in_name,strlen(in_name)+1); + + this->mask = xNew(Min); + xMemCpy(this->mask, maskin, Min); + + this->M=Min; + + } + /*}}}*/ + ~Regionaloutput(){/*{{{*/ + if(this->name)xDelete(this->name); + if(this->outputname)xDelete(this->outputname); + if(this->mask)xDelete(this->mask); + } + /*}}}*/ + /*Object virtual function resolutoin: */ + Object* copy() {/*{{{*/ + Regionaloutput* mf = new Regionaloutput(this->name,this->definitionenum,this->outputname,this->mask,this->M); + return (Object*) mf; + } + /*}}}*/ + void DeepEcho(void){/*{{{*/ + this->Echo(); + } + /*}}}*/ + void Echo(void){/*{{{*/ + _printf_(" Regionaloutput: " << this->name << " " << this->definitionenum << "\n"); + _printf_(" outputname enum: " << this->outputname << "Enum\n"); + _printf_(" mask: " << this->mask << "\n"); + _printf_(" M: " << this->M << "\n"); + } + /*}}}*/ + int Id(void){/*{{{*/ + return -1; + } + /*}}}*/ + void Marshall(char** pmarshalled_data,int* pmarshalled_data_size, int marshall_direction){/*{{{*/ + _error_("not implemented yet!"); + } + /*}}}*/ + int ObjectEnum(void){/*{{{*/ + return RegionaloutputEnum; + } + /*}}}*/ + /*Definition virtual function resolutoin: */ + int DefinitionEnum(){/*{{{*/ + + return this->definitionenum; + } + /*}}}*/ + char* Name(){/*{{{*/ + + char* name2=xNew(strlen(this->name)+1); + xMemCpy(name2,this->name,strlen(this->name)+1); + + return name2; + } + /*}}}*/ + IssmDouble Response(FemModel* femmodel){/*{{{*/ + + int i; + IssmDouble val_t=0.; + IssmDouble all_val_t=0.; + int outputenum = StringToEnumx(this->outputname); + + for(i=0;ielements->Size();i++){ + Element* element=(Element*)femmodel->elements->GetObjectByOffset(i); + switch(outputenum){ + case GroundedAreaEnum: + val_t+=element->GroundedArea(this->mask); + break; + case IceMassEnum: + val_t+=element->IceMass(this->mask); + break; + case IceVolumeEnum: + val_t+=element->IceVolume(this->mask); + break; + case IceVolumeAboveFloatationEnum: + val_t+=element->IceVolumeAboveFloatation(this->mask); + break; + case TotalFloatingBmbEnum: + val_t+=element->TotalFloatingBmb(this->mask); + break; + case TotalGroundedBmbEnum: + val_t+=element->TotalGroundedBmb(this->mask); + break; + case TotalSmbEnum: + val_t+=element->TotalSmb(this->mask); + break; + default: + _error_("Regional output type " << this->outputname << " not supported yet!"); + } + } + + ISSM_MPI_Allreduce ( (void*)&val_t,(void*)&all_val_t,1,ISSM_MPI_DOUBLE,ISSM_MPI_SUM,IssmComm::GetComm()); + val_t=all_val_t; + + return val_t; + } + /*}}}*/ +}; + +#endif /* _REGIONALOUTPUT_H_ */ Index: /issm/trunk-jpl/src/c/classes/classes.h =================================================================== --- /issm/trunk-jpl/src/c/classes/classes.h (revision 21807) +++ /issm/trunk-jpl/src/c/classes/classes.h (revision 21808) @@ -21,4 +21,5 @@ #include "./Masscon.h" #include "./Massconaxpby.h" +#include "./Regionaloutput.h" /*Constraints: */ Index: /issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateOutputDefinitions.cpp =================================================================== --- /issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateOutputDefinitions.cpp (revision 21807) +++ /issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateOutputDefinitions.cpp (revision 21808) @@ -183,18 +183,19 @@ /*masscon variables: */ int nummasscons; - char** masscon_name_s = NULL; - int* masscon_definitionenum_s = NULL; + char** masscon_name_s = NULL; + char** masscon_definitionstring_s = NULL; IssmDouble** masscon_levelset_s = NULL; - int* masscon_levelset_M_s = NULL; - int* masscon_levelset_N_s = NULL; + int* masscon_levelset_M_s = NULL; + int* masscon_levelset_N_s = NULL; /*Fetch name and levelset, etc ... (see src/m/classes/masscon.m): */ iomodel->FetchMultipleData(&masscon_name_s,&nummasscons, "md.masscon.name"); - iomodel->FetchMultipleData(&masscon_definitionenum_s,&nummasscons, "md.masscon.definitionenum"); + iomodel->FetchMultipleData(&masscon_definitionstring_s,&nummasscons, "md.masscon.definitionstring"); iomodel->FetchMultipleData(&masscon_levelset_s,&masscon_levelset_M_s,&masscon_levelset_N_s,&nummasscons,"md.masscon.levelset"); + for(j=0;jAddObject(new Masscon(masscon_name_s[j],masscon_definitionenum_s[j],masscon_levelset_s[j],masscon_levelset_M_s[j])); + output_definitions->AddObject(new Masscon(masscon_name_s[j],StringToEnumx(masscon_definitionstring_s[j]),masscon_levelset_s[j],masscon_levelset_M_s[j])); } @@ -206,4 +207,5 @@ string = masscon_name_s[j]; xDelete(string); + string = masscon_definitionstring_s[j]; xDelete(string); matrix = masscon_levelset_s[j]; xDelete(matrix); } @@ -212,5 +214,6 @@ xDelete(masscon_levelset_M_s); xDelete(masscon_levelset_N_s); - xDelete(masscon_definitionenum_s); + xDelete(masscon_definitionstring_s); + /*}}}*/ } @@ -257,4 +260,46 @@ /*}}}*/ } + else if (output_definition_enums[i]==RegionaloutputEnum){ + /*Deal with regional output: {{{*/ + + /*masscon variables: */ + int numout; + char** reg_name_s = NULL; + char** reg_definitionstring_s = NULL; + char** reg_outputnamestring_s = NULL; + IssmDouble** reg_mask_s = NULL; + int* reg_mask_M_s = NULL; + int* reg_mask_N_s = NULL; + + /*Fetch name and mask, etc ... (see src/m/classes/regionaloutput.m): */ + iomodel->FetchMultipleData(®_name_s,&numout, "md.regionaloutput.name"); + iomodel->FetchMultipleData(®_definitionstring_s,&numout, "md.regionaloutput.definitionstring"); + iomodel->FetchMultipleData(®_outputnamestring_s,&numout, "md.regionaloutput.outputnamestring"); + iomodel->FetchMultipleData(®_mask_s,®_mask_M_s,®_mask_N_s,&numout, "md.regionaloutput.mask"); + for(j=0;jAddObject(new Regionaloutput(reg_name_s[j],StringToEnumx(reg_definitionstring_s[j]),reg_outputnamestring_s[j],reg_mask_s[j],reg_mask_M_s[j])); + + } + + /*Free ressources:*/ + for(j=0;j(string); + string = reg_definitionstring_s[j]; xDelete(string); + string = reg_outputnamestring_s[j]; xDelete(string); + matrix = reg_mask_s[j]; xDelete(matrix); + } + xDelete(reg_name_s); + xDelete(reg_mask_s); + xDelete(reg_mask_M_s); + xDelete(reg_mask_N_s); + xDelete(reg_outputnamestring_s); + xDelete(reg_definitionstring_s); + } + /*}}}*/ else _error_("output definition enum " << output_definition_enums[i] << " not supported yet!"); } Index: /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h (revision 21807) +++ /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h (revision 21808) @@ -703,4 +703,5 @@ MassfluxatgateEnum, NodalvalueEnum, + RegionaloutputEnum, VxObsEnum, VyObsEnum, Index: /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp (revision 21807) +++ /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp (revision 21808) @@ -697,4 +697,5 @@ case MassfluxatgateEnum : return "Massfluxatgate"; case NodalvalueEnum : return "Nodalvalue"; + case RegionaloutputEnum : return "Regionaloutput"; case VxObsEnum : return "VxObs"; case VyObsEnum : return "VyObs"; Index: /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp (revision 21807) +++ /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp (revision 21808) @@ -712,4 +712,5 @@ else if (strcmp(name,"Massfluxatgate")==0) return MassfluxatgateEnum; else if (strcmp(name,"Nodalvalue")==0) return NodalvalueEnum; + else if (strcmp(name,"Regionaloutput")==0) return RegionaloutputEnum; else if (strcmp(name,"VxObs")==0) return VxObsEnum; else if (strcmp(name,"VyObs")==0) return VyObsEnum; Index: /issm/trunk-jpl/src/m/classes/masscon.m =================================================================== --- /issm/trunk-jpl/src/m/classes/masscon.m (revision 21807) +++ /issm/trunk-jpl/src/m/classes/masscon.m (revision 21808) @@ -3,10 +3,10 @@ % Usage: % masscon=masscon(); -% masscon=masscon('name','MassCon58',... %name of a North-East Greenland JPL MassCon +% masscon=masscon('name','MassCon58','definitionstring','Outputdefinition1',... %name of a North-East Greenland JPL MassCon % 'levelset',level); % % where level is a levelset vectorial field. % -% See also: MISFIT +% See also: MISFIT, MASSCONAXPBY, REGIONALOUTPUT classdef masscon @@ -14,9 +14,14 @@ %masscon name = ''; - definitionstring = ''; %string that identifies this output definition uniquely, from 'Outputdefinition[1-10]' + definitionstring = ''; %string that identifies this output definition uniquely, from 'Outputdefinition[1-10]' levelset = NaN; %levelset vectorial field which identifies the boundaries of the masscon end methods + function self = extrude(self,md) % {{{ + if ~isnan(self.levelset) + self.levelset=project3d(md,'vector',self.levelset,'type','node'); + end + end % }}} function self = masscon(varargin) % {{{ if nargin==0, @@ -40,5 +45,5 @@ error('masscon error message: ''name'' field should be a string!'); end - + OutputdefinitionStringArray={}; for i=1:100 @@ -46,11 +51,11 @@ end - md = checkfield(md,'fieldname','self.definitionstring','field',self.definitionstring,'values','Outputdefinition'); + md = checkfield(md,'fieldname','self.definitionstring','field',self.definitionstring,'values',OutputdefinitionStringArray); md = checkfield(md,'fieldname','self.levelset','field',self.levelset,'timeseries',1,'NaN',1,'Inf',1); end % }}} function md = disp(self) % {{{ - - disp(sprintf(' Misfit:\n')); + + disp(sprintf(' Masscon:\n')); fielddisplay(self,'name','identifier for this masscon response'); @@ -61,7 +66,7 @@ function md = marshall(self,prefix,md,fid) % {{{ - WriteData(fid,prefix,'object',self,'fieldname','name','format','String'); - WriteData(fid,prefix,'object',self,'fieldname','definitionstring','format','String'); - WriteData(fid,prefix,'object',self,'fieldname','levelset','format','DoubleMat','mattype',1); + WriteData(fid,prefix,'data',self.name,'name','md.masscon.name','format','String'); + WriteData(fid,prefix,'data',self.definitionstring,'name','md.masscon.definitionstring','format','String'); + WriteData(fid,prefix,'data',self.levelset,'name','md.masscon.levelset','format','DoubleMat','mattype',1); end % }}} Index: /issm/trunk-jpl/src/m/classes/massconaxpby.m =================================================================== --- /issm/trunk-jpl/src/m/classes/massconaxpby.m (revision 21807) +++ /issm/trunk-jpl/src/m/classes/massconaxpby.m (revision 21808) @@ -3,5 +3,5 @@ % Usage: % massconaxpby=massconaxpby(); -% massconaxpby=massconaxpby('name','MassCon58+35','namex','MassCon58','alpha',.5,'namey','MassCon35','beta',.5); +% massconaxpby=massconaxpby('name','MassCon58+35','namex','MassCon58','alpha',.5,'namey','MassCon35','beta',.5,'definitionstring','Outputdefinition1'); % % where name is the name of the massconaxpby object, namex is the name of the first masscon, namey the name of the second masscon and alpha,beta the @@ -9,18 +9,20 @@ % of two masscons. % -% See also: MASSCON +% See also: MASSCON, REGIONALOUTPUT classdef massconaxpby properties (SetAccess=public) %masscon axpby - name = ''; + name = ''; definitionstring = ''; %String that identifies this output definition uniquely, from 'Outputdefinition[1-10]' namex = ''; namey = ''; alpha = NaN; - beta = NaN; + beta = NaN; end methods + function self = extrude(self,md) % {{{ + end % }}} function self = massconaxpby(varargin) % {{{ if nargin==0, @@ -65,5 +67,5 @@ function md = disp(self) % {{{ - disp(sprintf(' Misfit:\n')); + disp(sprintf(' Massconaxpby:\n')); fielddisplay(self,'name','name'); @@ -77,10 +79,11 @@ function md = marshall(self,prefix,md,fid) % {{{ - WriteData(fid,prefix,'data',self.name,'name','md.massconaxpby.name','format','String'); - WriteData(fid,prefix,'data',self.definitionstring,'name','md.massconaxpby.definitionstring','format','String'); - WriteData(fid,prefix,'data',self.namex,'name','md.massconaxpby.namex','format','String'); - WriteData(fid,prefix,'data',self.namey,'name','md.massconaxpby.namey','format','String'); - WriteData(fid,prefix,'data',self.alpha,'name','md.massconaxpby.alpha','format','Double'); - WriteData(fid,prefix,'data',self.beta,'name','md.massconaxpby.beta','format','Double'); + WriteData(fid,prefix,'data',self.name,'name','md.massconaxpby.name','format','String'); + WriteData(fid,prefix,'data',self.definitionstring,'name','md.massconaxpby.definitionstring','format','String'); + WriteData(fid,prefix,'data',self.namex,'name','md.massconaxpby.namex','format','String'); + WriteData(fid,prefix,'data',self.namey,'name','md.massconaxpby.namey','format','String'); + WriteData(fid,prefix,'data',self.alpha,'name','md.massconaxpby.alpha','format','Double'); + WriteData(fid,prefix,'data',self.beta,'name','md.massconaxpby.beta','format','Double'); + end % }}} end Index: /issm/trunk-jpl/src/m/classes/massfluxatgate.m =================================================================== --- /issm/trunk-jpl/src/m/classes/massfluxatgate.m (revision 21807) +++ /issm/trunk-jpl/src/m/classes/massfluxatgate.m (revision 21808) @@ -16,4 +16,6 @@ end methods + function self = extrude(self,md) % {{{ + end % }}} function self = massfluxatgate(varargin) % {{{ if nargin==0, @@ -56,5 +58,5 @@ fielddisplay(self,'name','identifier for this massfluxatgate response'); fielddisplay(self,'profilename','name of file (shapefile or argus file) defining a profile (or gate)'); - fielddisplay(self,'definitionstring','string that identifies this output definition uniquely, from ''Outputdefinition[1-10]'''); + fielddisplay(self,'definitionstring','string that identifies this output definition uniquely, from ''Outputdefinition[1-100]'''); end % }}} Index: /issm/trunk-jpl/src/m/classes/massfluxatgate.py =================================================================== --- /issm/trunk-jpl/src/m/classes/massfluxatgate.py (revision 21807) +++ /issm/trunk-jpl/src/m/classes/massfluxatgate.py (revision 21808) @@ -35,8 +35,11 @@ string=" Massfluxatgate:" string="%s\n%s"%(string,fielddisplay(self,'name','identifier for this massfluxatgate response')) - string="%s\n%s"%(string,fielddisplay(self,'definitionstring','string that identifies this output definition uniquely, from Outputdefinition[1-10]')) + string="%s\n%s"%(string,fielddisplay(self,'definitionstring','string that identifies this output definition uniquely, from Outputdefinition[1-100]')) string="%s\n%s"%(string,fielddisplay(self,'profilename','name of file (shapefile or argus file) defining a profile (or gate)')) return string #}}} + def extrude(self,md): # {{{ + return self + #}}} def setdefaultparameters(self): # {{{ return self @@ -49,6 +52,11 @@ if not isinstance(self.profilename, basestring): raise RuntimeError("massfluxatgate error message: 'profilename' field should be a string!") - - md = checkfield(md,'field',self.definitionstring,'values',['Outputdefinition1','Outputdefinition2','Outputdefinition3','Outputdefinition4','Outputdefinition5','Outputdefinition6','Outputdefinition7','Outputdefinition8','Outputdefinition9','Outputdefinition10']) + + OutputdefinitionStringArray=[] + for i in range(1,100): + x='Outputdefinition'+str(i) + OutputdefinitionStringArray.append(x) + + md = checkfield(md,'field',self.definitionstring,'values',OutputdefinitionStringArray) #check the profilename points to a file!: Index: /issm/trunk-jpl/src/m/classes/model.m =================================================================== --- /issm/trunk-jpl/src/m/classes/model.m (revision 21807) +++ /issm/trunk-jpl/src/m/classes/model.m (revision 21808) @@ -260,6 +260,6 @@ if(md.hydrology.isefficientlayer==1) md.hydrology.spcepl_head=project2d(md,md.hydrology.spcepl_head,1); - end - end + end + end %materials @@ -308,4 +308,18 @@ if ~isnan(md.mask.ice_levelset), md.mask.ice_levelset=project2d(md,md.mask.ice_levelset,1); + end + + %outputdefinitions + for i=1:length(md.outputdefinition.definitions) + if isobject(md.outputdefinition.definitions{i}) + %get subfields + solutionsubfields=fields(md.outputdefinition.definitions{i}); + for j=1:length(solutionsubfields), + field=md.outputdefinition.definitions{i}.(solutionsubfields{j}); + if length(field)==md.mesh.numberofvertices | length(field)==md.mesh.numberofelements, + md.outputdefinition.definitions{i}.(solutionsubfields{j})=project2d(md,md.outputdefinition.definitions{i}.(solutionsubfields{j}),1); + end + end + end end @@ -607,4 +621,20 @@ end + %OutputDefinitions fields + for i=1:length(md1.outputdefinition.definitions), + if isobject(md1.outputdefinition.definitions{i}) + %get subfields + solutionsubfields=fields(md1.outputdefinition.definitions{i}); + for j=1:length(solutionsubfields), + field=md1.outputdefinition.definitions{i}.(solutionsubfields{j}); + if length(field)==numberofvertices1, + md2.outputdefinition.definitions{i}.(solutionsubfields{j})=field(pos_node); + elseif length(field)==numberofelements1, + md2.outputdefinition.definitions{i}.(solutionsubfields{j})=field(pos_elem); + end + end + end + end + %Keep track of pos_node and pos_elem md2.mesh.extractedvertices=pos_node; @@ -793,4 +823,5 @@ md.qmu=extrude(md.qmu,md); md.basalforcings=extrude(md.basalforcings,md); + md.outputdefinition=extrude(md.outputdefinition,md); %increase connectivity if less than 25: Index: /issm/trunk-jpl/src/m/classes/model.py =================================================================== --- /issm/trunk-jpl/src/m/classes/model.py (revision 21807) +++ /issm/trunk-jpl/src/m/classes/model.py (revision 21808) @@ -467,4 +467,21 @@ setattr(fieldr,solutionsubfield,subfield) + #OutputDefinitions fields + if md1.outputdefinition.definitions: + for solutionfield,field in md1.outputdefinition.__dict__.iteritems(): + if isinstance(field,list): + #get each definition + for i,fieldi in enumerate(field): + if fieldi: + fieldr=getattr(md2.outputdefinition,solutionfield)[i] + #get subfields + for solutionsubfield,subfield in fieldi.__dict__.iteritems(): + if np.size(subfield)==numberofvertices1: + setattr(fieldr,solutionsubfield,subfield[pos_node]) + elif np.size(subfield)==numberofelements1: + setattr(fieldr,solutionsubfield,subfield[pos_elem]) + else: + setattr(fieldr,solutionsubfield,subfield) + #Keep track of pos_node and pos_elem md2.mesh.extractedvertices=pos_node+1 @@ -649,4 +666,5 @@ md.qmu.extrude(md) md.basalforcings.extrude(md) + md.outputdefinition.extrude(md) #increase connectivity if less than 25: @@ -714,12 +732,12 @@ - # Hydrologydc variables - if hasattr(md.hydrology,'hydrologydc'): - md.hydrology.spcsediment_head=project2d(md,md.hydrology.spcsediment_head,1) - md.hydrology.mask_eplactive_node=project2d(md,md.hydrology.mask_eplactive_node,1) - md.hydrology.sediment_transmitivity=project2d(md,md.hydrology.sediment_transmitivity,1) - md.hydrology.basal_moulin_input=project2d(md,md.hydrology.basal_moulin_input,1) - if md.hydrology.isefficientlayer == 1: - md.hydrology.spcepl_head=project2d(md,md.hydrology.spcepl_head,1) + # Hydrologydc variables + if hasattr(md.hydrology,'hydrologydc'): + md.hydrology.spcsediment_head=project2d(md,md.hydrology.spcsediment_head,1) + md.hydrology.mask_eplactive_node=project2d(md,md.hydrology.mask_eplactive_node,1) + md.hydrology.sediment_transmitivity=project2d(md,md.hydrology.sediment_transmitivity,1) + md.hydrology.basal_moulin_input=project2d(md,md.hydrology.basal_moulin_input,1) + if md.hydrology.isefficientlayer == 1: + md.hydrology.spcepl_head=project2d(md,md.hydrology.spcepl_head,1) #boundary conditions @@ -761,4 +779,21 @@ md.mask.groundedice_levelset=project2d(md,md.mask.groundedice_levelset,1) md.mask.ice_levelset=project2d(md,md.mask.ice_levelset,1) + + #OutputDefinitions + if md.outputdefinition.definitions: + for solutionfield,field in md.outputdefinition.__dict__.iteritems(): + if isinstance(field,list): + #get each definition + for i,fieldi in enumerate(field): + if fieldi: + fieldr=getattr(md.outputdefinition,solutionfield)[i] + #get subfields + for solutionsubfield,subfield in fieldi.__dict__.iteritems(): + if np.size(subfield)==md.mesh.numberofvertices: + setattr(fieldr,solutionsubfield,project2d(md,subfield,1)) + elif np.size(subfield)==md.mesh.numberofelements: + setattr(fieldr,solutionsubfield,project2d(md,subfield,1)) + else: + setattr(fieldr,solutionsubfield,subfield) #Initialize with the 2d mesh Index: /issm/trunk-jpl/src/m/classes/outputdefinition.m =================================================================== --- /issm/trunk-jpl/src/m/classes/outputdefinition.m (revision 21807) +++ /issm/trunk-jpl/src/m/classes/outputdefinition.m (revision 21808) @@ -9,4 +9,9 @@ end methods + function self = extrude(self,md) % {{{ + for i=1:length(self.definitions) + self.definitions{i}=extrude(self.definitions{i},md); + end + end % }}} function self = outputdefinition(varargin) % {{{ switch nargin Index: /issm/trunk-jpl/src/m/classes/outputdefinition.py =================================================================== --- /issm/trunk-jpl/src/m/classes/outputdefinition.py (revision 21807) +++ /issm/trunk-jpl/src/m/classes/outputdefinition.py (revision 21808) @@ -22,4 +22,10 @@ return string #}}} + def extrude(self,md): # {{{ + for definition in self.definitions: + definition.extrude(md); + + return self + #}}} def setdefaultparameters(self): # {{{ return self Index: /issm/trunk-jpl/src/m/classes/regionaloutput.m =================================================================== --- /issm/trunk-jpl/src/m/classes/regionaloutput.m (revision 21808) +++ /issm/trunk-jpl/src/m/classes/regionaloutput.m (revision 21808) @@ -0,0 +1,106 @@ +%REGIONALOUTPUT class definition +% +% Usage: +% regionaloutput=regionaloutput(); +% regionaloutput=regionaloutput('name','Volume1','definitionstring','Outputdefinition1','outputnamestring','IceVolume',... +% 'mask',mask); +% regionaloutput=regionaloutput('name','Volume1','definitionstring','Outputdefinition1','outputnamestring','IceVolume',... +% 'maskexpstring','Exp/Mask.exp','model',md) +% +% where mask is a vectorial field of size md.mesh.numberofvertices,1 : where vertices with values > 1 are to be included in the calculated region. +% Alternatively, the user can pass in an Argus file and model object instead of a mask, and mask will be calculated for the user +% +% See also: MISFIT, MASSCON, MASSCONAXPBY + +classdef regionaloutput + properties (SetAccess=public) + %regionaloutput + name = ''; + definitionstring = ''; %string that identifies this output definition uniquely, from 'Outputdefinition[1-100]' + outputnamestring = ''; %string that identifies the type of output you want, eg. IceVolume, TotalSmb, GroudedArea + mask = NaN; %mask vectorial field which identifies the region of interest (value > 0 will be included) + maskexpstring = ''; + end + + methods + function self = extrude(self,md) % {{{ + if ~isnan(self.mask) + self.mask=project3d(md,'vector',self.mask,'type','node'); + end + end % }}} + function self = regionaloutput(varargin) % {{{ + if nargin==0, + self=setdefaultparameters(self); + else + %use provided options to change fields + options=pairoptions(varargin{:}); + + %get name + self.name=getfieldvalue(options,'name',''); + if nargin==8 + self.mask=getfieldvalue(options,'mask',NaN); + if isnan(self.mask) + error('regionaloutput error message: ''mask'' field or ''maskexpstring'' and ''model'' fields should be defined!'); + end + elseif nargin==10 + modelname=getfieldvalue(options,'model'); + self.maskexpstring=getfieldvalue(options,'maskexpstring'); + self=setmaskfromexp(self,modelname); + else + error('regionaloutput error message: ''mask'' field or ''maskexpstring'' and ''model'' fields should be defined!'); + end + + self.definitionstring=getfieldvalue(options,'definitionstring'); + self.outputnamestring=getfieldvalue(options,'outputnamestring'); + + end + end % }}} + function self = setdefaultparameters(self) % {{{ + end % }}} + function self = setmaskfromexp(self,md) % {{{ + + if length(self.maskexpstring)>0 + self.mask=ContourToMesh(md.mesh.elements,md.mesh.x,md.mesh.y,self.maskexpstring,'node',1); + end + end % }}} + function md = checkconsistency(self,md,solution,analyses) % {{{ + + if ~ischar(self.name), + error('regionaloutput error message: ''name'' field should be a string!'); + end + if ~ischar(self.outputnamestring), + error('regionaloutput error message: ''outputnamestring'' field should be a string!'); + end + + OutputdefinitionStringArray={}; + for i=1:100 + OutputdefinitionStringArray{i}=strcat('Outputdefinition',num2str(i)); + end + self=setmaskfromexp(self,md); + + md = checkfield(md,'fieldname','self.definitionstring','field',self.definitionstring,'values',OutputdefinitionStringArray); + md = checkfield(md,'fieldname','self.mask','field',self.mask,'size',[md.mesh.numberofvertices 1],'NaN',1,'Inf',1); + + end % }}} + function md = disp(self) % {{{ + + disp(sprintf(' Regionaloutput:\n')); + + fielddisplay(self,'name','identifier for this regional response'); + fielddisplay(self,'definitionstring','string that identifies this output definition uniquely, from ''Outputdefinition[1-100]'''); + fielddisplay(self,'outputnamestring','string that identifies the type of output you want, eg. IceVolume, TotalSmb, GroudedArea'); + fielddisplay(self,'mask','mask vectorial field which identifies the region of interest (value > 0 will be included)'); + fielddisplay(self,'maskexpstring','name of Argus file that can be passed in to define the regional mask'); + + end % }}} + function md = marshall(self,prefix,md,fid) % {{{ + + self=setmaskfromexp(self,md); + WriteData(fid,prefix,'data',self.name,'name','md.regionaloutput.name','format','String'); + WriteData(fid,prefix,'data',self.definitionstring,'name','md.regionaloutput.definitionstring','format','String'); + WriteData(fid,prefix,'data',self.outputnamestring,'name','md.regionaloutput.outputnamestring','format','String'); + WriteData(fid,prefix,'data',self.mask,'name','md.regionaloutput.mask','format','DoubleMat','mattype',1); + + end % }}} + end +end Index: /issm/trunk-jpl/src/m/classes/regionaloutput.py =================================================================== --- /issm/trunk-jpl/src/m/classes/regionaloutput.py (revision 21808) +++ /issm/trunk-jpl/src/m/classes/regionaloutput.py (revision 21808) @@ -0,0 +1,103 @@ +from fielddisplay import fielddisplay +from pairoptions import pairoptions +from checkfield import checkfield +from WriteData import WriteData +from MeshProfileIntersection import MeshProfileIntersection +import os + +class regionaloutput(object): + """ + REGIONALOUTPUT class definition + + Usage: + regionaloutput=regionaloutput(); + regionaloutput=regionaloutput('name','Volume1','definitionstring','Outputdefinition1','outputnamestring','IceVolume','mask',mask); + regionaloutput=regionaloutput('name','Volume1','definitionstring','Outputdefinition1','outputnamestring','IceVolume','maskexpstring','Exp/Mask.exp','model',md) + + where mask is a vectorial field of size md.mesh.numberofvertices,1 : where vertices with values > 1 are to be included in the calculated region. + Alternatively, the user can pass in an Argus file and model object instead of a mask, and mask will be calculated for the user + """ + + def __init__(self,*args): # {{{ + + self.name = '' + self.definitionstring = '' + self.outputnamestring = '' + self.mask = float('NaN') + self.maskexpstring = '' + + #set defaults + self.setdefaultparameters() + + #use provided options to change fields + options=pairoptions(*args) + + #OK get other fields + self=options.AssignObjectFields(self) + + #get name + modelname=getfieldvalue(options,'model'); + self.maskexpstring=getfieldvalue(options,'maskexpstring'); + self.setmaskfromexp(self,md) + error('regionaloutput error message: ''mask'' field or ''maskexpstring'' and ''model'' fields should be defined!'); + + #}}} + def __repr__(self): # {{{ + + string=" Regionaloutput:" + string="%s\n%s"%(string,fielddisplay(self,'name','identifier for this regional response')) + string="%s\n%s"%(string,fielddisplay(self,'definitionstring','string that identifies this output definition uniquely, from Outputdefinition[1-100]')) + string="%s\n%s"%(string,fielddisplay(self,'outputnamestring','string that identifies the type of output you want, eg. IceVolume, TotalSmb, GroudedArea')) + string="%s\n%s"%(string,fielddisplay(self,'mask','mask vectorial field which identifies the region of interest (value > 0 will be included)')) + string="%s\n%s"%(string,fielddisplay(self,'maskexpstring','name of Argus file that can be passed in to define the regional mask')) + return string + #}}} + def extrude(self,md): # {{{ + self.mask=project3d(md,'vector',self.mask,'type','node') + return self + #}}} + def setdefaultparameters(self): # {{{ + return self + #}}} + def setmaskfromexp(self,md): # {{{ + if len(self.maskexpstring,0) > 0: + self.mask=ContourToMesh(md.mesh.elements,md.mesh.x,md.mesh.y,self.maskexpstring,'node',1) + + return self + # }}} + def checkconsistency(self,md,solution,analyses): # {{{ + + if not isinstance(self.name, basestring): + raise RuntimeError("regionaloutput error message: 'name' field should be a string!") + + if not isinstance(self.outputnamestring, basestring): + raise RuntimeError("regionaloutput error message: 'outputnamestring' field should be a string!") + + if len(self.maskexpstring,0) > 0: + if not os.path.isfile(self.profilename): + raise RuntimeError("regionaloutput error message: file name for mask exp does not point to a legitimate file on disk!") + else: + self.setmaskfromexp(self,md) + + OutputdefinitionStringArray=[] + for i in range(1,100): + x='Outputdefinition'+str(i) + OutputdefinitionStringArray.append(x) + + md = checkfield(md,'field',self.definitionstring,'values',OutputdefinitionStringArray]) + md = checkfield(md,'field',self.mask,'size',[md.mesh.numberofvertices 1],'NaN',1,'Inf',1) + + return md + # }}} + def marshall(self,prefix,md,fid): # {{{ + + #before marshalling, make sure mask is set: + self.setmaskfromexp(self,md) + + #ok, marshall strings and mask: + WriteData(fid,prefix,'data',self.name,'name','md.regionaloutput.name','format','String') + WriteData(fid,prefix,'data',self.definitionstring,'name','md.regionaloutput.definitionstring','format','String') + WriteData(fid,prefix,'data',self.outputnamestring,'name','md.regionaloutput.outputnamestring','format','String'); + WriteData(fid,prefix,'data',self.mask,'name','md.regionaloutput.mask','format','DoubleMat','mattype',1); + + # }}}