Index: /issm/trunk/src/c/EnumDefinitions/EnumDefinitions.h =================================================================== --- /issm/trunk/src/c/EnumDefinitions/EnumDefinitions.h (revision 5280) +++ /issm/trunk/src/c/EnumDefinitions/EnumDefinitions.h (revision 5281) @@ -170,4 +170,6 @@ BedSlopeYEnum, BoundaryEnum, + CmResponseEnum, + CmResponsesEnum, CmMaxDmpSlopeEnum, CmMaxDmpValueEnum, @@ -190,5 +192,4 @@ EpsVelEnum, FillEnum, - FitEnum, FractionIncrementEnum, FrictionEnum, @@ -197,4 +198,9 @@ KflagEnum, MassFluxEnum, + SurfaceAbsVelMisfitEnum, + SurfaceRelVelMisfitEnum, + SurfaceLogVelMisfitEnum, + SurfaceLogVxVyMisfitEnum, + SurfaceAverageVelMisfitEnum, MaxPenetrationEnum, MeanVelEnum, Index: /issm/trunk/src/c/EnumDefinitions/EnumToString.cpp =================================================================== --- /issm/trunk/src/c/EnumDefinitions/EnumToString.cpp (revision 5280) +++ /issm/trunk/src/c/EnumDefinitions/EnumToString.cpp (revision 5281) @@ -144,4 +144,6 @@ case BedSlopeYEnum : return "BedSlopeY"; case BoundaryEnum : return "Boundary"; + case CmResponseEnum : return "CmResponse"; + case CmResponsesEnum : return "CmResponses"; case CmMaxDmpSlopeEnum : return "CmMaxDmpSlope"; case CmMaxDmpValueEnum : return "CmMaxDmpValue"; @@ -164,5 +166,4 @@ case EpsVelEnum : return "EpsVel"; case FillEnum : return "Fill"; - case FitEnum : return "Fit"; case FractionIncrementEnum : return "FractionIncrement"; case FrictionEnum : return "Friction"; @@ -171,4 +172,9 @@ case KflagEnum : return "Kflag"; case MassFluxEnum : return "MassFlux"; + case SurfaceAbsVelMisfitEnum : return "SurfaceAbsVelMisfit"; + case SurfaceRelVelMisfitEnum : return "SurfaceRelVelMisfit"; + case SurfaceLogVelMisfitEnum : return "SurfaceLogVelMisfit"; + case SurfaceLogVxVyMisfitEnum : return "SurfaceLogVxVyMisfit"; + case SurfaceAverageVelMisfitEnum : return "SurfaceAverageVelMisfit"; case MaxPenetrationEnum : return "MaxPenetration"; case MeanVelEnum : return "MeanVel"; Index: /issm/trunk/src/c/EnumDefinitions/StringToEnum.cpp =================================================================== --- /issm/trunk/src/c/EnumDefinitions/StringToEnum.cpp (revision 5280) +++ /issm/trunk/src/c/EnumDefinitions/StringToEnum.cpp (revision 5281) @@ -142,4 +142,6 @@ else if (strcmp(name,"BedSlopeY")==0) return BedSlopeYEnum; else if (strcmp(name,"Boundary")==0) return BoundaryEnum; + else if (strcmp(name,"CmResponse")==0) return CmResponseEnum; + else if (strcmp(name,"CmResponses")==0) return CmResponsesEnum; else if (strcmp(name,"CmMaxDmpSlope")==0) return CmMaxDmpSlopeEnum; else if (strcmp(name,"CmMaxDmpValue")==0) return CmMaxDmpValueEnum; @@ -162,5 +164,4 @@ else if (strcmp(name,"EpsVel")==0) return EpsVelEnum; else if (strcmp(name,"Fill")==0) return FillEnum; - else if (strcmp(name,"Fit")==0) return FitEnum; else if (strcmp(name,"FractionIncrement")==0) return FractionIncrementEnum; else if (strcmp(name,"Friction")==0) return FrictionEnum; @@ -169,4 +170,9 @@ else if (strcmp(name,"Kflag")==0) return KflagEnum; else if (strcmp(name,"MassFlux")==0) return MassFluxEnum; + else if (strcmp(name,"SurfaceAbsVelMisfit")==0) return SurfaceAbsVelMisfitEnum; + else if (strcmp(name,"SurfaceRelVelMisfit")==0) return SurfaceRelVelMisfitEnum; + else if (strcmp(name,"SurfaceLogVelMisfit")==0) return SurfaceLogVelMisfitEnum; + else if (strcmp(name,"SurfaceLogVxVyMisfit")==0) return SurfaceLogVxVyMisfitEnum; + else if (strcmp(name,"SurfaceAverageVelMisfit")==0) return SurfaceAverageVelMisfitEnum; else if (strcmp(name,"MaxPenetration")==0) return MaxPenetrationEnum; else if (strcmp(name,"MeanVel")==0) return MeanVelEnum; @@ -305,5 +311,2 @@ else if (strcmp(name,"WaitOnLock")==0) return WaitOnLockEnum; else if (strcmp(name,"Yts")==0) return YtsEnum; - else ISSMERROR("Enum %s not found",name); - -} Index: /issm/trunk/src/c/Makefile.am =================================================================== --- /issm/trunk/src/c/Makefile.am (revision 5280) +++ /issm/trunk/src/c/Makefile.am (revision 5281) @@ -433,6 +433,14 @@ ./modules/SurfaceAreax/SurfaceAreax.h\ ./modules/SurfaceAreax/SurfaceAreax.cpp\ - ./modules/Misfitx/Misfitx.h\ - ./modules/Misfitx/Misfitx.cpp\ + ./modules/SurfaceAbsVelMisfitx/SurfaceAbsVelMisfitx.h\ + ./modules/SurfaceAbsVelMisfitx/SurfaceAbsVelMisfitx.cpp\ + ./modules/SurfaceRelVelMisfitx/SurfaceRelVelMisfitx.h\ + ./modules/SurfaceRelVelMisfitx/SurfaceRelVelMisfitx.cpp\ + ./modules/SurfaceLogVelMisfitx/SurfaceLogVelMisfitx.h\ + ./modules/SurfaceLogVelMisfitx/SurfaceLogVelMisfitx.cpp\ + ./modules/SurfaceLogVxVyMisfitx/SurfaceLogVxVyMisfitx.h\ + ./modules/SurfaceLogVxVyMisfitx/SurfaceLogVxVyMisfitx.cpp\ + ./modules/SurfaceAverageVelMisfitx/SurfaceAverageVelMisfitx.h\ + ./modules/SurfaceAverageVelMisfitx/SurfaceAverageVelMisfitx.cpp\ ./modules/CostFunctionx/CostFunctionx.h\ ./modules/CostFunctionx/CostFunctionx.cpp\ @@ -968,6 +976,14 @@ ./modules/SurfaceAreax/SurfaceAreax.h\ ./modules/SurfaceAreax/SurfaceAreax.cpp\ - ./modules/Misfitx/Misfitx.h\ - ./modules/Misfitx/Misfitx.cpp\ + ./modules/SurfaceAbsVelMisfitx/SurfaceAbsVelMisfitx.h\ + ./modules/SurfaceAbsVelMisfitx/SurfaceAbsVelMisfitx.cpp\ + ./modules/SurfaceRelVelMisfitx/SurfaceRelVelMisfitx.h\ + ./modules/SurfaceRelVelMisfitx/SurfaceRelVelMisfitx.cpp\ + ./modules/SurfaceLogVelMisfitx/SurfaceLogVelMisfitx.h\ + ./modules/SurfaceLogVelMisfitx/SurfaceLogVelMisfitx.cpp\ + ./modules/SurfaceLogVxVyMisfitx/SurfaceLogVxVyMisfitx.h\ + ./modules/SurfaceLogVxVyMisfitx/SurfaceLogVxVyMisfitx.cpp\ + ./modules/SurfaceAverageVelMisfitx/SurfaceAverageVelMisfitx.h\ + ./modules/SurfaceAverageVelMisfitx/SurfaceAverageVelMisfitx.cpp\ ./modules/CostFunctionx/CostFunctionx.h\ ./modules/CostFunctionx/CostFunctionx.cpp\ Index: /issm/trunk/src/c/modules/CostFunctionx/CostFunctionx.cpp =================================================================== --- /issm/trunk/src/c/modules/CostFunctionx/CostFunctionx.cpp (revision 5280) +++ /issm/trunk/src/c/modules/CostFunctionx/CostFunctionx.cpp (revision 5281) @@ -9,7 +9,7 @@ #include "../../toolkits/toolkits.h" #include "../../EnumDefinitions/EnumDefinitions.h" -#include "../SurfaceAreax/SurfaceAreax.h" +#include "../Responsex/Responsex.h" -void CostFunctionx( double* pJ, Elements* elements,Nodes* nodes, Vertices* vertices, Loads* loads, Materials* materials,Parameters* parameters,bool process_units){ +void CostFunctionx( double* pJ, Elements* elements,Nodes* nodes, Vertices* vertices, Loads* loads, Materials* materials,Parameters* parameters,int response){ /*Intermediary*/ @@ -20,26 +20,21 @@ /*output: */ - double J=0; - double J_sum; + double J; + double Jreg=0; + double Jreg_sum; + + /*Get response*/ + Responsex(&J,elements,nodes,vertices,loads,materials,parameters,(const char*)EnumToString(response),false); //False means DO NOT process units - /*Compute surface area: */ - SurfaceAreax(&S,elements,nodes,vertices, loads,materials,parameters); - - /*add surface area to elements :*/ - for(i=0;iSize();i++){ - Element* element=(Element*)elements->GetObjectByOffset(i); - element->InputUpdateFromVector(&S,SurfaceAreaEnum,ConstantEnum); - } - - /*Compute cost function: */ + /*Add Regularization terms: */ for (i=0;iSize();i++){ element=(Element*)elements->GetObjectByOffset(i); - J+=element->CostFunction(process_units); + Jreg+=element->RegularizeInversion(); } /*Sum all J from all cpus of the cluster:*/ - MPI_Reduce (&J,&J_sum,1,MPI_DOUBLE,MPI_SUM,0,MPI_COMM_WORLD ); - MPI_Bcast(&J_sum,1,MPI_DOUBLE,0,MPI_COMM_WORLD); - J=J_sum; + MPI_Reduce (&Jreg,&Jreg_sum,1,MPI_DOUBLE,MPI_SUM,0,MPI_COMM_WORLD ); + MPI_Bcast(&Jreg_sum,1,MPI_DOUBLE,0,MPI_COMM_WORLD); + J=J+Jreg_sum; /*Assign output pointers: */ Index: /issm/trunk/src/c/modules/CostFunctionx/CostFunctionx.h =================================================================== --- /issm/trunk/src/c/modules/CostFunctionx/CostFunctionx.h (revision 5280) +++ /issm/trunk/src/c/modules/CostFunctionx/CostFunctionx.h (revision 5281) @@ -10,6 +10,5 @@ /* local prototypes: */ -void CostFunctionx( double* pJ, Elements* elements,Nodes* nodes, Vertices* vertices, Loads* loads, Materials* materials, Parameters* parameters,bool process_units); +void CostFunctionx( double* pJ, Elements* elements,Nodes* nodes, Vertices* vertices, Loads* loads, Materials* materials, Parameters* parameters,int response); -#endif /* _MISFITX_H */ - +#endif Index: /issm/trunk/src/c/modules/MassFluxx/MassFluxx.cpp =================================================================== --- /issm/trunk/src/c/modules/MassFluxx/MassFluxx.cpp (revision 5280) +++ /issm/trunk/src/c/modules/MassFluxx/MassFluxx.cpp (revision 5281) @@ -10,5 +10,5 @@ #include "../../EnumDefinitions/EnumDefinitions.h" -void MassFluxx(double* pmass_flux, Elements* elements,Nodes* nodes, Vertices* vertices,Loads* loads,Materials* materials, Parameters* parameters,char* descriptor,bool process_units){ +void MassFluxx(double* pmass_flux, Elements* elements,Nodes* nodes, Vertices* vertices,Loads* loads,Materials* materials, Parameters* parameters,const char* descriptor,bool process_units){ int i,j; Index: /issm/trunk/src/c/modules/MassFluxx/MassFluxx.h =================================================================== --- /issm/trunk/src/c/modules/MassFluxx/MassFluxx.h (revision 5280) +++ /issm/trunk/src/c/modules/MassFluxx/MassFluxx.h (revision 5281) @@ -10,5 +10,5 @@ /* local prototypes: */ -void MassFluxx(double* pmass_flux, Elements* elements,Nodes* nodes, Vertices* vertices,Loads* loads,Materials* materials, Parameters* parameters,char* descriptor,bool process_units); +void MassFluxx(double* pmass_flux, Elements* elements,Nodes* nodes, Vertices* vertices,Loads* loads,Materials* materials, Parameters* parameters,const char* descriptor,bool process_units); Index: /issm/trunk/src/c/modules/ModelProcessorx/Control/CreateParametersControl.cpp =================================================================== --- /issm/trunk/src/c/modules/ModelProcessorx/Control/CreateParametersControl.cpp (revision 5280) +++ /issm/trunk/src/c/modules/ModelProcessorx/Control/CreateParametersControl.cpp (revision 5281) @@ -60,15 +60,15 @@ /*Now, recover fit, optscal and maxiter as vectors: */ - IoModelFetchData(&iomodel->fit,NULL,NULL,iomodel_handle,"fit"); + IoModelFetchData(&iomodel->cm_responses,NULL,NULL,iomodel_handle,"cm_responses"); IoModelFetchData(&iomodel->cm_jump,NULL,NULL,iomodel_handle,"cm_jump"); IoModelFetchData(&iomodel->optscal,NULL,NULL,iomodel_handle,"optscal"); IoModelFetchData(&iomodel->maxiter,NULL,NULL,iomodel_handle,"maxiter"); - parameters->AddObject(new DoubleVecParam(FitEnum,iomodel->fit,iomodel->nsteps)); + parameters->AddObject(new DoubleVecParam(CmResponsesEnum,iomodel->cm_responses,iomodel->nsteps)); parameters->AddObject(new DoubleVecParam(CmJumpEnum,iomodel->cm_jump,iomodel->nsteps)); parameters->AddObject(new DoubleVecParam(OptScalEnum,iomodel->optscal,iomodel->nsteps)); parameters->AddObject(new DoubleVecParam(MaxIterEnum,iomodel->maxiter,iomodel->nsteps)); - xfree((void**)&iomodel->fit); + xfree((void**)&iomodel->cm_responses); xfree((void**)&iomodel->cm_jump); xfree((void**)&iomodel->optscal); Index: /issm/trunk/src/c/modules/ModelProcessorx/Qmu/CreateParametersQmu.cpp =================================================================== --- /issm/trunk/src/c/modules/ModelProcessorx/Qmu/CreateParametersQmu.cpp (revision 5280) +++ /issm/trunk/src/c/modules/ModelProcessorx/Qmu/CreateParametersQmu.cpp (revision 5281) @@ -207,13 +207,4 @@ qmu_mass_flux_present=true; } - - if (strcmp(descriptor,"misfit")==0){ - - /*We need the fit: */ - IoModelFetchData(&iomodel->fit,NULL,NULL,iomodel_handle,"fit"); - parameters->SetParam(iomodel->fit,iomodel->nsteps,FitEnum); - xfree((void**)&iomodel->fit); - - } } Index: /issm/trunk/src/c/modules/Responsex/Responsex.cpp =================================================================== --- /issm/trunk/src/c/modules/Responsex/Responsex.cpp (revision 5280) +++ /issm/trunk/src/c/modules/Responsex/Responsex.cpp (revision 5281) @@ -16,5 +16,5 @@ #include "../modules.h" -void Responsex(double* presponse,Elements* elements,Nodes* nodes, Vertices* vertices,Loads* loads,Materials* materials, Parameters* parameters,char* response_descriptor,bool process_units){ +void Responsex(double* presponse,Elements* elements,Nodes* nodes, Vertices* vertices,Loads* loads,Materials* materials, Parameters* parameters,const char* response_descriptor,bool process_units){ if(strcmp(response_descriptor,"MinVel")==0){ @@ -51,12 +51,24 @@ MaxAbsVzx( presponse, elements,nodes, vertices, loads, materials, parameters,process_units); } - else if(strcmp(response_descriptor,"Misfit")==0){ - Misfitx( presponse, elements,nodes, vertices, loads, materials, parameters,process_units); + else if(strcmp(response_descriptor,"SurfaceAbsVelMisfit")==0){ + SurfaceAbsVelMisfitx( presponse, elements,nodes, vertices, loads, materials, parameters,process_units); } - else if(strlen(response_descriptor)>=8){ - if(strncmp(response_descriptor,"MassFlux",8)==0) MassFluxx(presponse,elements,nodes,vertices,loads,materials,parameters,response_descriptor,process_units); + else if(strcmp(response_descriptor,"SurfaceRelVelMisfit")==0){ + SurfaceRelVelMisfitx( presponse, elements,nodes, vertices, loads, materials, parameters,process_units); + } + else if(strcmp(response_descriptor,"SurfaceLogVelMisfit")==0){ + SurfaceLogVelMisfitx( presponse, elements,nodes, vertices, loads, materials, parameters,process_units); + } + else if(strcmp(response_descriptor,"SurfaceLogVxVyMisfit")==0){ + SurfaceLogVxVyMisfitx( presponse, elements,nodes, vertices, loads, materials, parameters,process_units); + } + else if(strcmp(response_descriptor,"SurfaceAverageVelMisfit")==0){ + SurfaceAverageVelMisfitx( presponse, elements,nodes, vertices, loads, materials, parameters,process_units); + } + else if(strncmp(response_descriptor,"MassFlux",8)==0){ + MassFluxx(presponse,elements,nodes,vertices,loads,materials,parameters,response_descriptor,process_units); } else{ - ISSMERROR("%s%s%s"," response descriptor : ",response_descriptor," not supported yet!"); + ISSMERROR(" response descriptor \"%s\" not supported yet!",response_descriptor); } } Index: /issm/trunk/src/c/modules/Responsex/Responsex.h =================================================================== --- /issm/trunk/src/c/modules/Responsex/Responsex.h (revision 5280) +++ /issm/trunk/src/c/modules/Responsex/Responsex.h (revision 5281) @@ -9,5 +9,5 @@ #include "../../Container/Container.h" -void Responsex(double* presponse,Elements* elements,Nodes* nodes, Vertices* vertices,Loads* loads,Materials* materials, Parameters* parameters,char* response_descriptor,bool process_units); +void Responsex(double* presponse,Elements* elements,Nodes* nodes, Vertices* vertices,Loads* loads,Materials* materials, Parameters* parameters,const char* response_descriptor,bool process_units); #endif /* _RESPONSESXX_H */ Index: /issm/trunk/src/c/modules/modules.h =================================================================== --- /issm/trunk/src/c/modules/modules.h (revision 5280) +++ /issm/trunk/src/c/modules/modules.h (revision 5281) @@ -55,5 +55,9 @@ #include "./MinVyx/MinVyx.h" #include "./MinVzx/MinVzx.h" -#include "./Misfitx/Misfitx.h" +#include "./SurfaceAbsVelMisfitx/SurfaceAbsVelMisfitx.h" +#include "./SurfaceRelVelMisfitx/SurfaceRelVelMisfitx.h" +#include "./SurfaceLogVelMisfitx/SurfaceLogVelMisfitx.h" +#include "./SurfaceLogVxVyMisfitx/SurfaceLogVxVyMisfitx.h" +#include "./SurfaceAverageVelMisfitx/SurfaceAverageVelMisfitx.h" #include "./ModelProcessorx/ModelProcessorx.h" #include "./NodeConnectivityx/NodeConnectivityx.h" Index: /issm/trunk/src/c/objects/Elements/Element.h =================================================================== --- /issm/trunk/src/c/objects/Elements/Element.h (revision 5280) +++ /issm/trunk/src/c/objects/Elements/Element.h (revision 5281) @@ -42,6 +42,10 @@ virtual void GradjDrag(Vec gradient)=0; virtual void GradjB(Vec gradient)=0; - virtual double Misfit(bool process_units)=0; - virtual double CostFunction(bool process_units)=0; + virtual double SurfaceAbsVelMisfit(bool process_units)=0; + virtual double SurfaceRelVelMisfit(bool process_units)=0; + virtual double SurfaceLogVelMisfit(bool process_units)=0; + virtual double SurfaceLogVxVyMisfit(bool process_units)=0; + virtual double SurfaceAverageVelMisfit(bool process_units)=0; + virtual double RegularizeInversion(void)=0; virtual double SurfaceArea(void)=0; virtual void InputDepthAverageAtBase(int enum_type,int average_enum_type,int object_enum)=0; Index: /issm/trunk/src/c/objects/Elements/Penta.cpp =================================================================== --- /issm/trunk/src/c/objects/Elements/Penta.cpp (revision 5280) +++ /issm/trunk/src/c/objects/Elements/Penta.cpp (revision 5281) @@ -666,6 +666,6 @@ } /*}}}*/ -/*FUNCTION Penta::CostFunction {{{1*/ -double Penta::CostFunction(bool process_units){ +/*FUNCTION Penta::RegularizeInversion {{{1*/ +double Penta::RegularizeInversion(void){ double J; @@ -696,5 +696,5 @@ /*This element should be collapsed into a tria element at its base. Create this tria element, - * and compute CostFunction*/ + * and compute RegularizeInversion*/ /*Depth Average B*/ @@ -702,5 +702,5 @@ tria=(Tria*)SpawnTria(0,1,2); //grids 0, 1 and 2 make the new tria (lower face). - J=tria->CostFunction(process_units); + J=tria->RegularizeInversion(); delete tria->matice; delete tria; @@ -720,5 +720,5 @@ tria=(Tria*)SpawnTria(3,4,5); //grids 3, 4 and 5 make the new tria (upper face). - J=tria->CostFunction(process_units); + J=tria->RegularizeInversion(); delete tria->matice; delete tria; @@ -1597,6 +1597,6 @@ } /*}}}*/ -/*FUNCTION Penta::Misfit {{{1*/ -double Penta::Misfit(bool process_units){ +/*FUNCTION Penta::SurfaceAbsVelMisfit {{{1*/ +double Penta::SurfaceAbsVelMisfit(bool process_units){ double J; @@ -1627,7 +1627,7 @@ /*This element should be collapsed into a tria element at its base. Create this tria element, - * and compute Misfit*/ + * and compute SurfaceAbsVelMisfit*/ tria=(Tria*)SpawnTria(0,1,2); //grids 0, 1 and 2 make the new tria (lower face). - J=tria->Misfit(process_units); + J=tria->SurfaceAbsVelMisfit(process_units); delete tria->matice; delete tria; return J; @@ -1636,5 +1636,185 @@ tria=(Tria*)SpawnTria(3,4,5); //grids 3, 4 and 5 make the new tria (upper face). - J=tria->Misfit(process_units); + J=tria->SurfaceAbsVelMisfit(process_units); + delete tria->matice; delete tria; + return J; + } +} +/*}}}*/ +/*FUNCTION Penta::SurfaceRelVelMisfit {{{1*/ +double Penta::SurfaceRelVelMisfit(bool process_units){ + + double J; + Tria* tria=NULL; + + /*inputs: */ + bool onwater; + bool onsurface; + bool onbed; + int approximation; + + /*retrieve inputs :*/ + inputs->GetParameterValue(&onwater,ElementOnWaterEnum); + inputs->GetParameterValue(&onbed,ElementOnBedEnum); + inputs->GetParameterValue(&onsurface,ElementOnSurfaceEnum); + inputs->GetParameterValue(&approximation,ApproximationEnum); + + /*If on water, return 0: */ + if(onwater)return 0; + + /*Bail out if this element if: + * -> Non MacAyeal and not on the surface + * -> MacAyeal (2d model) and not on bed) */ + if ((approximation!=MacAyealApproximationEnum && !onsurface) || (MacAyealApproximationEnum && !onbed)){ + return 0; + } + else if (MacAyealApproximationEnum){ + + /*This element should be collapsed into a tria element at its base. Create this tria element, + * and compute SurfaceRelVelMisfit*/ + tria=(Tria*)SpawnTria(0,1,2); //grids 0, 1 and 2 make the new tria (lower face). + J=tria->SurfaceRelVelMisfit(process_units); + delete tria->matice; delete tria; + return J; + } + else{ + + tria=(Tria*)SpawnTria(3,4,5); //grids 3, 4 and 5 make the new tria (upper face). + J=tria->SurfaceRelVelMisfit(process_units); + delete tria->matice; delete tria; + return J; + } +} +/*}}}*/ +/*FUNCTION Penta::SurfaceLogVelMisfit {{{1*/ +double Penta::SurfaceLogVelMisfit(bool process_units){ + + double J; + Tria* tria=NULL; + + /*inputs: */ + bool onwater; + bool onsurface; + bool onbed; + int approximation; + + /*retrieve inputs :*/ + inputs->GetParameterValue(&onwater,ElementOnWaterEnum); + inputs->GetParameterValue(&onbed,ElementOnBedEnum); + inputs->GetParameterValue(&onsurface,ElementOnSurfaceEnum); + inputs->GetParameterValue(&approximation,ApproximationEnum); + + /*If on water, return 0: */ + if(onwater)return 0; + + /*Bail out if this element if: + * -> Non MacAyeal and not on the surface + * -> MacAyeal (2d model) and not on bed) */ + if ((approximation!=MacAyealApproximationEnum && !onsurface) || (MacAyealApproximationEnum && !onbed)){ + return 0; + } + else if (MacAyealApproximationEnum){ + + /*This element should be collapsed into a tria element at its base. Create this tria element, + * and compute SurfaceLogVelMisfit*/ + tria=(Tria*)SpawnTria(0,1,2); //grids 0, 1 and 2 make the new tria (lower face). + J=tria->SurfaceLogVelMisfit(process_units); + delete tria->matice; delete tria; + return J; + } + else{ + + tria=(Tria*)SpawnTria(3,4,5); //grids 3, 4 and 5 make the new tria (upper face). + J=tria->SurfaceLogVelMisfit(process_units); + delete tria->matice; delete tria; + return J; + } +} +/*}}}*/ +/*FUNCTION Penta::SurfaceLogVxVyMisfit {{{1*/ +double Penta::SurfaceLogVxVyMisfit(bool process_units){ + + double J; + Tria* tria=NULL; + + /*inputs: */ + bool onwater; + bool onsurface; + bool onbed; + int approximation; + + /*retrieve inputs :*/ + inputs->GetParameterValue(&onwater,ElementOnWaterEnum); + inputs->GetParameterValue(&onbed,ElementOnBedEnum); + inputs->GetParameterValue(&onsurface,ElementOnSurfaceEnum); + inputs->GetParameterValue(&approximation,ApproximationEnum); + + /*If on water, return 0: */ + if(onwater)return 0; + + /*Bail out if this element if: + * -> Non MacAyeal and not on the surface + * -> MacAyeal (2d model) and not on bed) */ + if ((approximation!=MacAyealApproximationEnum && !onsurface) || (MacAyealApproximationEnum && !onbed)){ + return 0; + } + else if (MacAyealApproximationEnum){ + + /*This element should be collapsed into a tria element at its base. Create this tria element, + * and compute SurfaceLogVxVyMisfit*/ + tria=(Tria*)SpawnTria(0,1,2); //grids 0, 1 and 2 make the new tria (lower face). + J=tria->SurfaceLogVxVyMisfit(process_units); + delete tria->matice; delete tria; + return J; + } + else{ + + tria=(Tria*)SpawnTria(3,4,5); //grids 3, 4 and 5 make the new tria (upper face). + J=tria->SurfaceLogVxVyMisfit(process_units); + delete tria->matice; delete tria; + return J; + } +} +/*}}}*/ +/*FUNCTION Penta::SurfaceAverageVelMisfit {{{1*/ +double Penta::SurfaceAverageVelMisfit(bool process_units){ + + double J; + Tria* tria=NULL; + + /*inputs: */ + bool onwater; + bool onsurface; + bool onbed; + int approximation; + + /*retrieve inputs :*/ + inputs->GetParameterValue(&onwater,ElementOnWaterEnum); + inputs->GetParameterValue(&onbed,ElementOnBedEnum); + inputs->GetParameterValue(&onsurface,ElementOnSurfaceEnum); + inputs->GetParameterValue(&approximation,ApproximationEnum); + + /*If on water, return 0: */ + if(onwater)return 0; + + /*Bail out if this element if: + * -> Non MacAyeal and not on the surface + * -> MacAyeal (2d model) and not on bed) */ + if ((approximation!=MacAyealApproximationEnum && !onsurface) || (MacAyealApproximationEnum && !onbed)){ + return 0; + } + else if (MacAyealApproximationEnum){ + + /*This element should be collapsed into a tria element at its base. Create this tria element, + * and compute SurfaceAverageVelMisfit*/ + tria=(Tria*)SpawnTria(0,1,2); //grids 0, 1 and 2 make the new tria (lower face). + J=tria->SurfaceAverageVelMisfit(process_units); + delete tria->matice; delete tria; + return J; + } + else{ + + tria=(Tria*)SpawnTria(3,4,5); //grids 3, 4 and 5 make the new tria (upper face). + J=tria->SurfaceAverageVelMisfit(process_units); delete tria->matice; delete tria; return J; @@ -5900,5 +6080,5 @@ name==TemperatureEnum || name==ControlParameterEnum || - name==FitEnum || + name==CmResponseEnum || name==DragCoefficientEnum || name==GradientEnum || Index: /issm/trunk/src/c/objects/Elements/Penta.h =================================================================== --- /issm/trunk/src/c/objects/Elements/Penta.h (revision 5280) +++ /issm/trunk/src/c/objects/Elements/Penta.h (revision 5281) @@ -71,5 +71,5 @@ void Configure(Elements* elements,Loads* loads,DataSet* nodes,Materials* materials,Parameters* parameters); void SetCurrentConfiguration(Elements* elements,Loads* loads,DataSet* nodes,Materials* materials,Parameters* parameters); - double CostFunction(bool process_units); + double RegularizeInversion(void); void CreateKMatrix(Mat Kgg); void CreatePVector(Vec pg); @@ -104,5 +104,9 @@ void MinVy(double* pminvy, bool process_units); void MinVz(double* pminvz, bool process_units); - double Misfit(bool process_units); + double SurfaceAbsVelMisfit(bool process_units); + double SurfaceRelVelMisfit(bool process_units); + double SurfaceLogVelMisfit(bool process_units); + double SurfaceLogVxVyMisfit(bool process_units); + double SurfaceAverageVelMisfit(bool process_units); void PatchFill(int* pcount, Patch* patch); void PatchSize(int* pnumrows, int* pnumvertices,int* pnumnodes); Index: /issm/trunk/src/c/objects/Elements/Tria.cpp =================================================================== --- /issm/trunk/src/c/objects/Elements/Tria.cpp (revision 5280) +++ /issm/trunk/src/c/objects/Elements/Tria.cpp (revision 5281) @@ -548,11 +548,11 @@ } /*}}}*/ -/*FUNCTION Tria::CostFunction {{{1*/ -double Tria::CostFunction(bool process_units){ +/*FUNCTION Tria::RegularizeInversion {{{1*/ +double Tria::RegularizeInversion(){ int i; /* output: */ - double Jelem; + double Jelem=0; /* node data: */ @@ -602,7 +602,4 @@ /* Get node coordinates and dof list: */ GetVerticesCoordinates(&xyz_list[0][0], nodes, numgrids); - - /*First, get Misfit*/ - Jelem=Misfit(process_units); /* Get gaussian points and weights (make this a statically initialized list of points? fstd): */ @@ -1860,6 +1857,6 @@ } /*}}}*/ -/*FUNCTION Tria::Misfit {{{1*/ -double Tria::Misfit(bool process_units){ +/*FUNCTION Tria::SurfaceAbsVelMisfit {{{1*/ +double Tria::SurfaceAbsVelMisfit(bool process_units){ int i; @@ -1900,10 +1897,4 @@ double Jdet; - /*relative and logarithmic control method :*/ - double scalex=1; - double scaley=1; - double S=0; - int fit=-1; - /*inputs: */ bool onwater; @@ -1921,8 +1912,4 @@ /* Recover input data: */ - inputs->GetParameterValue(&fit,FitEnum); - if(fit==3){ - inputs->GetParameterValue(&S,SurfaceAreaEnum); - } inputs->GetParameterValues(&obs_vx_list[0],&gaussgrids[0][0],3,VxObsEnum); inputs->GetParameterValues(&obs_vy_list[0],&gaussgrids[0][0],3,VyObsEnum); @@ -1935,5 +1922,5 @@ this->parameters->FindParam(&epsvel,EpsVelEnum); - /* Compute Misfit at the 3 nodes + /* Compute SurfaceAbsVelMisfit at the 3 nodes * Here we integrate linearized functions: * @@ -1944,94 +1931,17 @@ * the vertex i */ - if(fit==0){ - /*We are using an absolute misfit: - * - * 1 [ 2 2 ] - * J = --- | (u - u ) + (v - v ) | - * 2 [ obs obs ] - * - */ - for (i=0;iparameters); - - } - else if(fit==1){ - /*We are using a relative misfit: - * - * 1 [ \bar{v}^2 2 \bar{v}^2 2 ] - * J = --- | ------------- (u - u ) + ------------- (v - v ) | - * 2 [ (u + eps)^2 obs (v + eps)^2 obs ] - * obs obs - */ - for (i=0;iparameters); - - } - else if(fit==2){ - /*We are using a logarithmic misfit: - * - * [ vel + eps ] 2 - * J = 4 \bar{v}^2 | log ( ----------- ) | - * [ vel + eps ] - * obs - */ - for (i=0;iparameters); - - } - else if(fit==3){ - /*We are using a spacially average absolute misfit: - * - * 1 2 2 - * J = --- sqrt( (u - u ) + (v - v ) ) - * S obs obs - */ - for (i=0;iparameters); - - /*Take the square root, and scale by surface: */ - for (i=0;iparameters); - - } - else{ - /*Not supported yet! : */ - ISSMERROR("%s%g","unsupported type of fit: ",fit); - } + + /*We are using an absolute misfit: + * + * 1 [ 2 2 ] + * J = --- | (u - u ) + (v - v ) | + * 2 [ obs obs ] + * + */ + for (i=0;iparameters); /*Apply weights to misfits*/ @@ -2057,5 +1967,5 @@ TriaRef::GetParameterValue(&misfit, &misfit_list[0],gauss_l1l2l3); - /*compute Misfit*/ + /*compute SurfaceAbsVelMisfit*/ Jelem+=misfit*Jdet*gauss_weight; } @@ -2066,4 +1976,529 @@ xfree((void**)&gauss_weights); + /*Return: */ + return Jelem; +} +/*}}}*/ +/*FUNCTION Tria::SurfaceRelVelMisfit {{{1*/ +double Tria::SurfaceRelVelMisfit(bool process_units){ + + int i; + + /* output: */ + double Jelem=0; + + /* node data: */ + const int numgrids=3; + const int numdof=2*numgrids; + const int NDOF2=2; + double xyz_list[numgrids][3]; + + /* grid data: */ + double vx_list[numgrids]; + double vy_list[numgrids]; + double obs_vx_list[numgrids]; + double obs_vy_list[numgrids]; + double misfit_square_list[numgrids]; + double misfit_list[numgrids]; + double weights_list[numgrids]; + + /* gaussian points: */ + int num_gauss,ig; + double* first_gauss_area_coord = NULL; + double* second_gauss_area_coord = NULL; + double* third_gauss_area_coord = NULL; + double* gauss_weights = NULL; + double gauss_weight; + double gauss_l1l2l3[3]; + double gaussgrids[numgrids][numgrids]={{1,0,0},{0,1,0},{0,0,1}}; + + /* parameters: */ + double velocity_mag,obs_velocity_mag; + double misfit; + + /* Jacobian: */ + double Jdet; + + /*relative and logarithmic control method :*/ + double scalex=1; + double scaley=1; + + /*inputs: */ + bool onwater; + + double meanvel, epsvel; + + /*retrieve inputs :*/ + inputs->GetParameterValue(&onwater,ElementOnWaterEnum); + + /*If on water, return 0: */ + if(onwater)return 0; + + /* Get node coordinates and dof list: */ + GetVerticesCoordinates(&xyz_list[0][0], nodes, numgrids); + + /* Recover input data: */ + inputs->GetParameterValues(&obs_vx_list[0],&gaussgrids[0][0],3,VxObsEnum); + inputs->GetParameterValues(&obs_vy_list[0],&gaussgrids[0][0],3,VyObsEnum); + inputs->GetParameterValues(&vx_list[0],&gaussgrids[0][0],3,VxEnum); + inputs->GetParameterValues(&vy_list[0],&gaussgrids[0][0],3,VyEnum); + inputs->GetParameterValues(&weights_list[0],&gaussgrids[0][0],3,WeightsEnum); + + /*retrieve some parameters: */ + this->parameters->FindParam(&meanvel,MeanVelEnum); + this->parameters->FindParam(&epsvel,EpsVelEnum); + + /* Compute SurfaceRelVelMisfit at the 3 nodes + * Here we integrate linearized functions: + * + * J(E) = int_E sum_{i=1}^3 J_i Phi_i + * + * where J_i are the misfits at the 3 nodes of the triangle + * Phi_i is the nodal function (P1) with respect to + * the vertex i + */ + + /*We are using a relative misfit: + * + * 1 [ \bar{v}^2 2 \bar{v}^2 2 ] + * J = --- | ------------- (u - u ) + ------------- (v - v ) | + * 2 [ (u + eps)^2 obs (v + eps)^2 obs ] + * obs obs + */ + for (i=0;iparameters); + + /*Apply weights to misfits*/ + for (i=0;iGetParameterValue(&onwater,ElementOnWaterEnum); + + /*If on water, return 0: */ + if(onwater)return 0; + + /* Get node coordinates and dof list: */ + GetVerticesCoordinates(&xyz_list[0][0], nodes, numgrids); + + /* Recover input data: */ + inputs->GetParameterValues(&obs_vx_list[0],&gaussgrids[0][0],3,VxObsEnum); + inputs->GetParameterValues(&obs_vy_list[0],&gaussgrids[0][0],3,VyObsEnum); + inputs->GetParameterValues(&vx_list[0],&gaussgrids[0][0],3,VxEnum); + inputs->GetParameterValues(&vy_list[0],&gaussgrids[0][0],3,VyEnum); + inputs->GetParameterValues(&weights_list[0],&gaussgrids[0][0],3,WeightsEnum); + + /*retrieve some parameters: */ + this->parameters->FindParam(&meanvel,MeanVelEnum); + this->parameters->FindParam(&epsvel,EpsVelEnum); + + /* Compute SurfaceLogVelMisfit at the 3 nodes + * Here we integrate linearized functions: + * + * J(E) = int_E sum_{i=1}^3 J_i Phi_i + * + * where J_i are the misfits at the 3 nodes of the triangle + * Phi_i is the nodal function (P1) with respect to + * the vertex i + */ + + /*We are using a logarithmic misfit: + * + * [ vel + eps ] 2 + * J = 4 \bar{v}^2 | log ( ----------- ) | + * [ vel + eps ] + * obs + */ + for (i=0;iparameters); + + /*Apply weights to misfits*/ + for (i=0;iGetParameterValue(&onwater,ElementOnWaterEnum); + + /*If on water, return 0: */ + if(onwater)return 0; + + /* Get node coordinates and dof list: */ + GetVerticesCoordinates(&xyz_list[0][0], nodes, numgrids); + + /* Recover input data: */ + inputs->GetParameterValues(&obs_vx_list[0],&gaussgrids[0][0],3,VxObsEnum); + inputs->GetParameterValues(&obs_vy_list[0],&gaussgrids[0][0],3,VyObsEnum); + inputs->GetParameterValues(&vx_list[0],&gaussgrids[0][0],3,VxEnum); + inputs->GetParameterValues(&vy_list[0],&gaussgrids[0][0],3,VyEnum); + inputs->GetParameterValues(&weights_list[0],&gaussgrids[0][0],3,WeightsEnum); + + /*retrieve some parameters: */ + this->parameters->FindParam(&meanvel,MeanVelEnum); + this->parameters->FindParam(&epsvel,EpsVelEnum); + + /* Compute SurfaceLogVxVyMisfit at the 3 nodes + * Here we integrate linearized functions: + * + * J(E) = int_E sum_{i=1}^3 J_i Phi_i + * + * where J_i are the misfits at the 3 nodes of the triangle + * Phi_i is the nodal function (P1) with respect to + * the vertex i + */ + + /*We are using an logarithmic 2 misfit: + * + * 1 [ |u| + eps 2 |v| + eps 2 ] + * J = --- \bar{v}^2 | log ( ----------- ) + log ( ----------- ) | + * 2 [ |u |+ eps |v |+ eps ] + * obs obs + */ + for (i=0;iparameters); + + /*Apply weights to misfits*/ + for (i=0;iGetParameterValue(&onwater,ElementOnWaterEnum); + + /*If on water, return 0: */ + if(onwater)return 0; + + /* Get node coordinates and dof list: */ + GetVerticesCoordinates(&xyz_list[0][0], nodes, numgrids); + + /* Recover input data: */ + inputs->GetParameterValue(&S,SurfaceAreaEnum); + inputs->GetParameterValues(&obs_vx_list[0],&gaussgrids[0][0],3,VxObsEnum); + inputs->GetParameterValues(&obs_vy_list[0],&gaussgrids[0][0],3,VyObsEnum); + inputs->GetParameterValues(&vx_list[0],&gaussgrids[0][0],3,VxEnum); + inputs->GetParameterValues(&vy_list[0],&gaussgrids[0][0],3,VyEnum); + inputs->GetParameterValues(&weights_list[0],&gaussgrids[0][0],3,WeightsEnum); + + /*retrieve some parameters: */ + this->parameters->FindParam(&meanvel,MeanVelEnum); + this->parameters->FindParam(&epsvel,EpsVelEnum); + + /* Compute SurfaceAverageVelMisfit at the 3 nodes + * Here we integrate linearized functions: + * + * J(E) = int_E sum_{i=1}^3 J_i Phi_i + * + * where J_i are the misfits at the 3 nodes of the triangle + * Phi_i is the nodal function (P1) with respect to + * the vertex i + */ + + /*We are using a spacially average absolute misfit: + * + * 1 2 2 + * J = --- sqrt( (u - u ) + (v - v ) ) + * S obs obs + */ + for (i=0;iparameters); + + /*Take the square root, and scale by surface: */ + for (i=0;iGetParameterValue(&fit,FitEnum); inputs->GetParameterValue(&onwater,ElementOnWaterEnum); - - /*If fit!=3, do not compute surface: */ - if(fit!=3)return 0; /*If on water, return 0: */ @@ -4471,5 +4901,5 @@ double scale=0; double S=0; - int fit=-1; + int response; /*retrieve some parameters: */ @@ -4490,6 +4920,6 @@ inputs->GetParameterValues(&weights_list[0],&gaussgrids[0][0],3,WeightsEnum); - inputs->GetParameterValue(&fit,FitEnum); - if(fit==3){ + inputs->GetParameterValue(&response,CmResponseEnum); + if(response==SurfaceAverageVelMisfitEnum){ inputs->GetParameterValue(&S,SurfaceAreaEnum); } @@ -4508,5 +4938,5 @@ * the vertex i */ - if(fit==0){ + if(response==SurfaceAbsVelMisfitEnum){ /*We are using an absolute misfit: * @@ -4524,5 +4954,5 @@ } } - else if(fit==1){ + else if(response==SurfaceRelVelMisfitEnum){ /*We are using a relative misfit: * @@ -4546,5 +4976,5 @@ } } - else if(fit==2){ + else if(response==SurfaceLogVelMisfitEnum){ /*We are using a logarithmic misfit: * @@ -4567,5 +4997,5 @@ } } - else if(fit==3){ + else if(response==SurfaceAverageVelMisfitEnum){ /*We are using an spacially average absolute misfit: * @@ -4584,5 +5014,5 @@ } } - else if(fit==4){ + else if(response==SurfaceLogVxVyMisfitEnum){ /*We are using an logarithmic 2 misfit: * @@ -4604,5 +5034,5 @@ else{ /*Not supported yet! : */ - ISSMERROR("%s%g","unsupported type of fit: ",fit); + ISSMERROR("response %s not supported yet",EnumToString(response)); } @@ -4710,4 +5140,5 @@ double S=0; int fit=-1; + int response; /*retrieve some parameters: */ @@ -4728,6 +5159,6 @@ inputs->GetParameterValues(&weights_list[0],&gaussgrids[0][0],3,WeightsEnum); - inputs->GetParameterValue(&fit,FitEnum); - if(fit==3){ + inputs->GetParameterValue(&response,CmResponseEnum); + if(response==SurfaceAverageVelMisfitEnum){ inputs->GetParameterValue(&S,SurfaceAreaEnum); } @@ -4746,5 +5177,5 @@ * the vertex i */ - if(fit==0){ + if(response==SurfaceAbsVelMisfitEnum){ /*We are using an absolute misfit: * @@ -4762,5 +5193,5 @@ } } - else if(fit==1){ + else if(response==SurfaceRelVelMisfitEnum){ /*We are using a relative misfit: * @@ -4784,5 +5215,5 @@ } } - else if(fit==2){ + else if(response==SurfaceLogVelMisfitEnum){ /*We are using a logarithmic misfit: * @@ -4805,5 +5236,5 @@ } } - else if(fit==3){ + else if(response==SurfaceAverageVelMisfitEnum){ /*We are using an spacially average absolute misfit: * @@ -4822,5 +5253,5 @@ } } - else if(fit==4){ + else if(response==SurfaceLogVxVyMisfitEnum){ /*We are using an logarithmic 2 misfit: * @@ -4842,5 +5273,5 @@ else{ /*Not supported yet! : */ - ISSMERROR("%s%g","unsupported type of fit: ",fit); + ISSMERROR("response %s not supported yet",EnumToString(response)); } @@ -6298,5 +6729,5 @@ name==VxObsEnum || name==VyObsEnum || - name==FitEnum || + name==CmResponseEnum || name==DragCoefficientEnum || name==GradientEnum || Index: /issm/trunk/src/c/objects/Elements/Tria.h =================================================================== --- /issm/trunk/src/c/objects/Elements/Tria.h (revision 5280) +++ /issm/trunk/src/c/objects/Elements/Tria.h (revision 5281) @@ -68,5 +68,5 @@ void Configure(Elements* elements,Loads* loads,DataSet* nodes,Materials* materials,Parameters* parameters); void SetCurrentConfiguration(Elements* elements,Loads* loads,DataSet* nodes,Materials* materials,Parameters* parameters); - double CostFunction(bool process_units); + double RegularizeInversion(void); void CreateKMatrix(Mat Kgg); void CreatePVector(Vec pg); @@ -103,5 +103,9 @@ void MinVy(double* pminvy, bool process_units); void MinVz(double* pminvz, bool process_units); - double Misfit(bool process_units); + double SurfaceAbsVelMisfit(bool process_units); + double SurfaceRelVelMisfit(bool process_units); + double SurfaceLogVelMisfit(bool process_units); + double SurfaceLogVxVyMisfit(bool process_units); + double SurfaceAverageVelMisfit(bool process_units); void PatchFill(int* pcount, Patch* patch); void PatchSize(int* pnumrows, int* pnumvertices,int* pnumnodes); Index: /issm/trunk/src/c/objects/IoModel.cpp =================================================================== --- /issm/trunk/src/c/objects/IoModel.cpp (revision 5280) +++ /issm/trunk/src/c/objects/IoModel.cpp (revision 5281) @@ -86,5 +86,5 @@ xfree((void**)&this->rheology_B); xfree((void**)&this->rheology_n); - xfree((void**)&this->fit); + xfree((void**)&this->cm_responses); xfree((void**)&this->weights); xfree((void**)&this->cm_jump); @@ -302,5 +302,5 @@ /*!solution parameters: */ - this->fit=NULL; + this->cm_responses=NULL; this->weights=NULL; this->cm_jump=NULL; Index: /issm/trunk/src/c/objects/IoModel.h =================================================================== --- /issm/trunk/src/c/objects/IoModel.h (revision 5280) +++ /issm/trunk/src/c/objects/IoModel.h (revision 5281) @@ -129,5 +129,5 @@ /*solution parameters: */ - double* fit; + double* cm_responses; double* weights; double* cm_jump; Index: /issm/trunk/src/c/solutions/control_core.cpp =================================================================== --- /issm/trunk/src/c/solutions/control_core.cpp (revision 5280) +++ /issm/trunk/src/c/solutions/control_core.cpp (revision 5281) @@ -26,5 +26,5 @@ int solution_type; - double* fit=NULL; + double* responses=NULL; double* optscal=NULL; double* maxiter=NULL; @@ -46,5 +46,5 @@ femmodel->parameters->FindParam(&nsteps,NStepsEnum); femmodel->parameters->FindParam(&control_type,ControlTypeEnum); - femmodel->parameters->FindParam(&fit,NULL,FitEnum); + femmodel->parameters->FindParam(&responses,NULL,CmResponsesEnum); femmodel->parameters->FindParam(&optscal,NULL,OptScalEnum); femmodel->parameters->FindParam(&maxiter,NULL,MaxIterEnum); @@ -65,5 +65,5 @@ } - /*Initialize misfit: */ + /*Initialize responses: */ J=(double*)xmalloc(nsteps*sizeof(double)); @@ -75,6 +75,7 @@ for(n=0;nelements,femmodel->nodes, femmodel->vertices, femmodel->loads, femmodel->materials, femmodel->parameters,(int)fit[n],FitEnum); + InputUpdateFromConstantx(femmodel->elements,femmodel->nodes, femmodel->vertices, femmodel->loads, femmodel->materials, femmodel->parameters,(int)responses[n],CmResponseEnum); /*In case we are running a steady state control method, compute new temperature field using new parameter distribution: */ @@ -97,5 +98,5 @@ InputControlUpdatex(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,search_scalar*optscal[n],true); - if(controlconvergence(J,fit,eps_cm,n)) break; + if(controlconvergence(J,responses,eps_cm,n)) break; /*Temporary saving every 5 control steps: */ @@ -117,5 +118,5 @@ cleanup_and_return: /*Free ressources: */ - xfree((void**)&fit); + xfree((void**)&responses); xfree((void**)&optscal); xfree((void**)&maxiter); Index: /issm/trunk/src/c/solutions/controlconvergence.cpp =================================================================== --- /issm/trunk/src/c/solutions/controlconvergence.cpp (revision 5280) +++ /issm/trunk/src/c/solutions/controlconvergence.cpp (revision 5281) @@ -18,5 +18,5 @@ -bool controlconvergence(double* J, double* fit, double eps_cm, int n){ +bool controlconvergence(double* J,double* responses, double eps_cm, int n){ int i; @@ -29,5 +29,5 @@ //go through the previous misfits(starting from n-2) while(i>=0){ - if (fit[i]==fit[n]){ + if (responses[i]==responses[n]){ //convergence test only if we have the same misfits if ((J[i]-J[n])/J[n] <= eps_cm){ Index: /issm/trunk/src/c/solutions/objectivefunctionC.cpp =================================================================== --- /issm/trunk/src/c/solutions/objectivefunctionC.cpp (revision 5280) +++ /issm/trunk/src/c/solutions/objectivefunctionC.cpp (revision 5281) @@ -28,15 +28,14 @@ /*parameters: */ - FemModel *femmodel = NULL; - int n; - double *optscal = NULL; - double *fit = NULL; - int control_type; - int solution_type; - int analysis_type; - bool isstokes = false; - bool conserve_loads = true; - bool process_units = false; - + FemModel *femmodel = NULL; + int n; + double *optscal = NULL; + double *responses =NULL; + int control_type; + int solution_type; + int analysis_type; + bool isstokes = false; + bool conserve_loads = true; + /*Recover finite element model: */ femmodel=optargs->femmodel; @@ -46,5 +45,5 @@ femmodel->parameters->FindParam(&optscal,NULL,OptScalEnum); - femmodel->parameters->FindParam(&fit,NULL,FitEnum); + femmodel->parameters->FindParam(&responses,NULL,CmResponsesEnum); femmodel->parameters->FindParam(&isstokes,IsStokesEnum); femmodel->parameters->FindParam(&control_type,ControlTypeEnum); @@ -82,9 +81,8 @@ /*Compute misfit for this velocity field.*/ - InputUpdateFromConstantx( femmodel->elements,femmodel->nodes, femmodel->vertices, femmodel->loads, femmodel->materials, femmodel->parameters,(int)fit[n],FitEnum); - CostFunctionx( &J, femmodel->elements,femmodel->nodes, femmodel->vertices,femmodel->loads, femmodel->materials, femmodel->parameters,process_units); + CostFunctionx( &J, femmodel->elements,femmodel->nodes, femmodel->vertices,femmodel->loads, femmodel->materials, femmodel->parameters,(int)responses[n]); /*Free ressources:*/ - xfree((void**)&fit); + xfree((void**)&responses); xfree((void**)&optscal); return J; Index: /issm/trunk/src/c/solutions/solutions.h =================================================================== --- /issm/trunk/src/c/solutions/solutions.h (revision 5280) +++ /issm/trunk/src/c/solutions/solutions.h (revision 5281) @@ -32,5 +32,5 @@ //convergence: void convergence(int* pconverged, Mat K_ff,Vec p_f,Vec u_f,Vec u_f_old,Parameters* parameters); -bool controlconvergence(double* J, double* fit, double eps_cm, int n); +bool controlconvergence(double* J, double* response, double eps_cm, int n); bool steadystateconvergence(FemModel* femmodel);