Index: /issm/trunk-jpl/src/c/classes/Misfit.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Misfit.h (revision 21737) +++ /issm/trunk-jpl/src/c/classes/Misfit.h (revision 21738) @@ -19,4 +19,6 @@ /*}}}*/ IssmDouble OutputDefinitionsResponsex(FemModel* femmodel,int output_enum); +void GetVectorFromInputsx( IssmDouble** pvector,FemModel* femmodel,int name,int type); + class Misfit: public Object, public Definition{ @@ -25,5 +27,5 @@ int definitionenum; - bool local; + int local; int model_enum; char* name; @@ -44,5 +46,5 @@ this->weights_enum = UNDEF; this->timeinterpolation=NULL; - this->local=true; + this->local=1; this->misfit=0; this->lock=0; @@ -50,5 +52,5 @@ } /*}}}*/ - Misfit(char* in_name, int in_definitionenum, int in_model_enum, int in_observation_enum, char* in_timeinterpolation, bool in_local, int in_weights_enum){/*{{{*/ + Misfit(char* in_name, int in_definitionenum, int in_model_enum, int in_observation_enum, char* in_timeinterpolation, int in_local, int in_weights_enum){/*{{{*/ this->definitionenum=in_definitionenum; @@ -133,5 +135,5 @@ femmodel->parameters->FindParam(&dt,TimesteppingTimeStepEnum); - if (this->local){ /*local computation: {{{*/ + if (this->local==1){ /*area integration using elements: {{{*/ int i; @@ -167,4 +169,70 @@ /*What we return is the value of misfit / time: */ return misfit/(time-starttime); + } /*}}}*/ + else if (this->local==2){ /*vertex by vertex computation: {{{*/ + + IssmDouble* model = NULL; + IssmDouble* observation= NULL; + IssmDouble* weights= NULL; + + /*Are we transient?:*/ + if (time==0){ + IssmDouble misfit_t=0.; + + /*get global vectors: */ + GetVectorFromInputsx(&model,femmodel,model_enum,VertexSIdEnum); + GetVectorFromInputsx(&observation,femmodel,observation_enum,VertexSIdEnum); + GetVectorFromInputsx(&weights,femmodel,weights_enum,VertexSIdEnum); + + int count=0; + for (int i=0;ivertices->NumberOfVertices();i++){ + misfit_t += pow(model[i]-observation[i],2)*weights[i]; + if (weights[i]!=0)count++; + } + misfit=sqrt(misfit_t)/count; + + /*Free ressources:*/ + xDelete(model); + xDelete(observation); + xDelete(weights); + + /*return value: */ + return misfit; + } + else{ + + IssmDouble misfit_t=0.; + IssmDouble all_misfit_t=0.; + + /*If we are locked, return time average: */ + if(this->lock)return misfit/(time-starttime); + + /*get global vectors: */ + GetVectorFromInputsx(&model,femmodel,model_enum,VertexSIdEnum); + GetVectorFromInputsx(&observation,femmodel,observation_enum,VertexSIdEnum); + GetVectorFromInputsx(&weights,femmodel,weights_enum,VertexSIdEnum); + + int count=0; + for (int i=0;ivertices->NumberOfVertices();i++){ + misfit_t += pow(model[i]-observation[i],2)*weights[i]; + if (weights[i]!=0)count++; + } + + /*Add this time's contribution to curent misfit: */ + misfit=sqrt(misfit_t)/count; + misfit+=dt*misfit_t; + + /*Do we lock? i.e. are we at final_time? :*/ + if(time==finaltime)this->lock=1; + + /*Free ressources:*/ + xDelete(model); + xDelete(observation); + xDelete(weights); + + /*What we return is the value of misfit / time: */ + return misfit/(time-starttime); + } + } /*}}}*/ else{ /*global computation: {{{ */ Index: /issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateOutputDefinitions.cpp =================================================================== --- /issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateOutputDefinitions.cpp (revision 21737) +++ /issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateOutputDefinitions.cpp (revision 21738) @@ -96,5 +96,5 @@ /*First create a misfit object for that specific string (misfit_model_string_s[j]):*/ - output_definitions->AddObject(new Misfit(misfit_name_s[j],StringToEnumx(misfit_definitionstring_s[j]),StringToEnumx(misfit_model_string_s[j]),StringToEnumx(misfit_observation_string_s[j]),misfit_timeinterpolation_s[j],(bool)misfit_local_s[j],StringToEnumx(misfit_weights_string_s[j]))); + output_definitions->AddObject(new Misfit(misfit_name_s[j],StringToEnumx(misfit_definitionstring_s[j]),StringToEnumx(misfit_model_string_s[j]),StringToEnumx(misfit_observation_string_s[j]),misfit_timeinterpolation_s[j],misfit_local_s[j],StringToEnumx(misfit_weights_string_s[j]))); /*Now, for this particular misfit object, make sure we plug into the elements: the observation, and the weights.*/ Index: /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h (revision 21737) +++ /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h (revision 21738) @@ -787,4 +787,6 @@ SealevelAbsoluteEnum, SealevelEustaticEnum, + SealevelObsEnum, + SealevelWeightsEnum, SealevelriseDeltathicknessEnum, SealevelriseMaxiterEnum, Index: /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp (revision 21737) +++ /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp (revision 21738) @@ -767,4 +767,6 @@ case SealevelAbsoluteEnum : return "SealevelAbsolute"; case SealevelEustaticEnum : return "SealevelEustatic"; + case SealevelObsEnum : return "SealevelObs"; + case SealevelWeightsEnum : return "SealevelWeights"; case SealevelriseDeltathicknessEnum : return "SealevelriseDeltathickness"; case SealevelriseMaxiterEnum : return "SealevelriseMaxiter"; Index: /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp (revision 21737) +++ /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp (revision 21738) @@ -785,4 +785,6 @@ else if (strcmp(name,"SealevelAbsolute")==0) return SealevelAbsoluteEnum; else if (strcmp(name,"SealevelEustatic")==0) return SealevelEustaticEnum; + else if (strcmp(name,"SealevelObs")==0) return SealevelObsEnum; + else if (strcmp(name,"SealevelWeights")==0) return SealevelWeightsEnum; else if (strcmp(name,"SealevelriseDeltathickness")==0) return SealevelriseDeltathicknessEnum; else if (strcmp(name,"SealevelriseMaxiter")==0) return SealevelriseMaxiterEnum; @@ -873,10 +875,10 @@ else if (strcmp(name,"SpcStatic")==0) return SpcStaticEnum; else if (strcmp(name,"SpcTransient")==0) return SpcTransientEnum; - else if (strcmp(name,"StringArrayParam")==0) return StringArrayParamEnum; - else if (strcmp(name,"StringParam")==0) return StringParamEnum; else stage=8; } if(stage==8){ - if (strcmp(name,"Seg")==0) return SegEnum; + if (strcmp(name,"StringArrayParam")==0) return StringArrayParamEnum; + else if (strcmp(name,"StringParam")==0) return StringParamEnum; + else if (strcmp(name,"Seg")==0) return SegEnum; else if (strcmp(name,"SegInput")==0) return SegInputEnum; else if (strcmp(name,"Tria")==0) return TriaEnum; @@ -996,10 +998,10 @@ else if (strcmp(name,"P2xP4")==0) return P2xP4Enum; else if (strcmp(name,"P1P1")==0) return P1P1Enum; - else if (strcmp(name,"P1P1GLS")==0) return P1P1GLSEnum; - else if (strcmp(name,"MINI")==0) return MINIEnum; else stage=9; } if(stage==9){ - if (strcmp(name,"MINIcondensed")==0) return MINIcondensedEnum; + if (strcmp(name,"P1P1GLS")==0) return P1P1GLSEnum; + else if (strcmp(name,"MINI")==0) return MINIEnum; + else if (strcmp(name,"MINIcondensed")==0) return MINIcondensedEnum; else if (strcmp(name,"TaylorHood")==0) return TaylorHoodEnum; else if (strcmp(name,"LATaylorHood")==0) return LATaylorHoodEnum; Index: /issm/trunk-jpl/src/m/classes/misfit.m =================================================================== --- /issm/trunk-jpl/src/m/classes/misfit.m (revision 21737) +++ /issm/trunk-jpl/src/m/classes/misfit.m (revision 21738) @@ -18,5 +18,5 @@ properties (SetAccess=public) %misfit - name = ''; + name = ''; definitionstring = ''; %string that identifies this output definition uniquely, from 'Outputdefinition[1-10]' model_string = ''; %string for field that is modeled @@ -24,5 +24,5 @@ observation_string = ''; %string for observed field. timeinterpolation = ''; - local = 1; + local = 1; weights = NaN; %weight coefficients for every vertex weights_string = ''; %string to identify this particular set of weights