Index: /issm/trunk-jpl/src/c/Makefile.am =================================================================== --- /issm/trunk-jpl/src/c/Makefile.am (revision 13698) +++ /issm/trunk-jpl/src/c/Makefile.am (revision 13699) @@ -318,8 +318,4 @@ ./modules/ConstraintsStatex/ConstraintsStatex.h\ ./modules/ConstraintsStatex/ConstraintsStateLocal.h\ - ./modules/Responsex/Responsex.h\ - ./modules/Responsex/Responsex.cpp\ - ./modules/RequestedOutputsx/RequestedOutputsx.h\ - ./modules/RequestedOutputsx/RequestedOutputsx.cpp\ ./modules/RequestedDependentsx/RequestedDependentsx.h\ ./modules/RequestedDependentsx/RequestedDependentsx.cpp\ @@ -363,6 +359,4 @@ ./modules/InputUpdateFromDakotax/InputUpdateFromDakotax.h\ ./modules/InputUpdateFromDakotax/InputUpdateFromDakotax.cpp\ - ./modules/DakotaResponsesx/DakotaResponsesx.h\ - ./modules/DakotaResponsesx/DakotaResponsesx.cpp\ ./modules/InputUpdateFromVectorDakotax/InputUpdateFromVectorDakotax.h\ ./modules/InputUpdateFromVectorDakotax/InputUpdateFromVectorDakotax.cpp\ @@ -437,6 +431,4 @@ ./modules/ThicknessAbsMisfitx/ThicknessAbsMisfitx.h\ ./modules/ThicknessAbsMisfitx/ThicknessAbsMisfitx.cpp\ - ./modules/CostFunctionx/CostFunctionx.h\ - ./modules/CostFunctionx/CostFunctionx.cpp\ ./modules/Orthx/Orthx.h\ ./modules/Orthx/Orthx.cpp\ @@ -445,6 +437,4 @@ ./modules/DragCoefficientAbsGradientx/DragCoefficientAbsGradientx.cpp\ ./modules/DragCoefficientAbsGradientx/DragCoefficientAbsGradientx.h\ - ./modules/ThicknessAbsGradientx/ThicknessAbsGradientx.cpp\ - ./modules/ThicknessAbsGradientx/ThicknessAbsGradientx.h\ ./modules/ThicknessAlongGradientx/ThicknessAlongGradientx.cpp\ ./modules/ThicknessAlongGradientx/ThicknessAlongGradientx.h\ Index: /issm/trunk-jpl/src/c/classes/FemModel.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/FemModel.cpp (revision 13698) +++ /issm/trunk-jpl/src/c/classes/FemModel.cpp (revision 13699) @@ -412,2 +412,260 @@ } /*}}}*/ +void FemModel::Responsex(IssmDouble* responses,const char* response_descriptor,bool process_units,int weight_index){/*{{{*/ + + + int response_descriptor_enum; + + response_descriptor_enum=StringToEnumx(response_descriptor); + this->Responsex(responses, response_descriptor_enum, process_units, weight_index); + +} +/*}}}*/ +void FemModel::Responsex(IssmDouble* responses,int response_descriptor_enum,bool process_units,int weight_index){/*{{{*/ + + + switch (response_descriptor_enum){ + + #ifdef _HAVE_RESPONSES_ + case IceVolumeEnum: IceVolumex( responses, elements,nodes, vertices, loads, materials, parameters,process_units); break; + case MinVelEnum: MinVelx( responses, elements,nodes, vertices, loads, materials, parameters,process_units); break; + case MaxVelEnum: MaxVelx( responses, elements,nodes, vertices, loads, materials, parameters,process_units); break; + case MinVxEnum: MinVxx( responses, elements,nodes, vertices, loads, materials, parameters,process_units); break; + case MaxVxEnum: MaxVxx( responses, elements,nodes, vertices, loads, materials, parameters,process_units); break; + case MaxAbsVxEnum: MaxAbsVxx( responses, elements,nodes, vertices, loads, materials, parameters,process_units); break; + case MinVyEnum: MinVyx( responses, elements,nodes, vertices, loads, materials, parameters,process_units); break; + case MaxVyEnum: MaxVyx( responses, elements,nodes, vertices, loads, materials, parameters,process_units); break; + case MaxAbsVyEnum: MaxAbsVyx( responses, elements,nodes, vertices, loads, materials, parameters,process_units); break; + case MinVzEnum: MinVzx( responses, elements,nodes, vertices, loads, materials, parameters,process_units); break; + case MaxVzEnum: MaxVzx( responses, elements,nodes, vertices, loads, materials, parameters,process_units); break; + case MaxAbsVzEnum: MaxAbsVzx( responses, elements,nodes, vertices, loads, materials, parameters,process_units); break; + case MassFluxEnum: MassFluxx( responses, elements,nodes, vertices, loads, materials, parameters,process_units); break; + case SurfaceAbsVelMisfitEnum: SurfaceAbsVelMisfitx( responses, elements,nodes, vertices, loads, materials, parameters,process_units,weight_index); break; + case SurfaceRelVelMisfitEnum: SurfaceRelVelMisfitx( responses, elements,nodes, vertices, loads, materials, parameters,process_units,weight_index); break; + case SurfaceLogVelMisfitEnum: SurfaceLogVelMisfitx( responses, elements,nodes, vertices, loads, materials, parameters,process_units,weight_index); break; + case SurfaceLogVxVyMisfitEnum: SurfaceLogVxVyMisfitx( responses, elements,nodes, vertices, loads, materials, parameters,process_units,weight_index); break; + case SurfaceAverageVelMisfitEnum:SurfaceAverageVelMisfitx( responses, elements,nodes, vertices, loads, materials, parameters,process_units,weight_index); break; + case ThicknessAbsMisfitEnum: ThicknessAbsMisfitx( responses, elements,nodes, vertices, loads, materials, parameters,process_units,weight_index); break; + case ThicknessAbsGradientEnum: this->ThicknessAbsGradientx( responses, process_units,weight_index); break; + case ThicknessAlongGradientEnum: ThicknessAlongGradientx( responses, elements,nodes, vertices, loads, materials, parameters,process_units,weight_index); break; + case ThicknessAcrossGradientEnum: ThicknessAcrossGradientx( responses, elements,nodes, vertices, loads, materials, parameters,process_units,weight_index); break; + case TotalSmbEnum: TotalSmbx( responses, elements,nodes, vertices, loads, materials, parameters,process_units); break; + case RheologyBbarAbsGradientEnum:RheologyBbarAbsGradientx( responses, elements,nodes, vertices, loads, materials, parameters,process_units,weight_index); break; + case DragCoefficientAbsGradientEnum:DragCoefficientAbsGradientx(responses, elements,nodes, vertices, loads, materials, parameters,process_units,weight_index); break; + case MaterialsRheologyBbarEnum:ElementResponsex(responses, elements,nodes, vertices, loads, materials, parameters,MaterialsRheologyBbarEnum,process_units); break; + case VelEnum:ElementResponsex(responses, elements,nodes, vertices, loads, materials, parameters,VelEnum,process_units); break; + case FrictionCoefficientEnum:NodalValuex(responses, FrictionCoefficientEnum,elements,nodes, vertices, loads, materials, parameters,process_units); break; + default: _error_("response descriptor \"" << EnumToStringx(response_descriptor_enum) << "\" not supported yet!"); break; + #else + default: _error_("ISSM was not compiled with responses capabilities, exiting!"); + #endif + } + +} +/*}}}*/ +void FemModel::RequestedOutputsx(int* requested_outputs, int numoutputs){/*{{{*/ + + + int output_enum; + int step; + IssmDouble time; + IssmDouble output_value; + Element *element = NULL; + + /*Get time and step*/ + parameters->FindParam(&step,StepEnum); + parameters->FindParam(&time,TimeEnum); + + /*retrieve Inputs*/ + if(numoutputs){ + for(int i=0;iAddObject(new GenericExternalResult(results->Size()+1,IceVolumeEnum,reCast(output_value),step,time)); + break; + case TotalSmbEnum: + Responsex(&output_value,"TotalSmb",false,0); + results->AddObject(new GenericExternalResult(results->Size()+1,TotalSmbEnum,reCast(output_value),step,time)); + break; + case MaxVelEnum: + Responsex(&output_value,"MaxVel",false,0); + results->AddObject(new GenericExternalResult(results->Size()+1,MaxVelEnum,reCast(output_value),step,time)); + break; + default: + /*create this output in the element inputs, and then transfer to results:*/ + for(int j=0;jSize();j++){ + element=(Element*)elements->GetObjectByOffset(j); + element->RequestedOutput(output_enum,step,time); + } + break; + } + } + } +} +/*}}}*/ +#ifdef _HAVE_CONTROL_ +void FemModel::ThicknessAbsGradientx( IssmDouble* pJ, bool process_units, int weight_index){/*{{{*/ + + + /*Intermediary*/ + int i; + Element* element=NULL; + + /*output: */ + IssmDouble J=0; + IssmDouble J_sum; + + /*Compute Misfit: */ + for (i=0;iSize();i++){ + element=(Element*)elements->GetObjectByOffset(i); + J+=element->ThicknessAbsGradient(process_units,weight_index); + } + + /*Sum all J from all cpus of the cluster:*/ + #ifdef _HAVE_MPI_ + MPI_Reduce (&J,&J_sum,1,MPI_DOUBLE,MPI_SUM,0,IssmComm::GetComm() ); + MPI_Bcast(&J_sum,1,MPI_DOUBLE,0,IssmComm::GetComm()); + J=J_sum; + #endif + + /*Assign output pointers: */ + *pJ=J; +} +/*}}}*/ +#endif + +#ifdef _HAVE_DAKOTA_ +void FemModel::CostFunctionx(IssmDouble* pJ){/*{{{*/ + + + /*Intermediary*/ + int i; + int num_responses; + Element *element = NULL; + int *responses = NULL; + + /*output: */ + IssmDouble J,Jplus; + + /*Recover parameters*/ + parameters->FindParam(&num_responses,InversionNumCostFunctionsEnum); + parameters->FindParam(&responses,NULL,NULL,StepResponsesEnum); + + /*Get response*/ + J=0; + for(int i=0;iResponsex(&Jplus,EnumToStringx(responses[i]),false,i); //False means DO NOT process units + J+=Jplus; + } + + /*Assign output pointers: */ + xDelete(responses); + *pJ=J; +} +/*}}}*/ +void FemModel::DakotaResponsesx(double* d_responses,char** responses_descriptors,int numresponsedescriptors,int d_numresponses){/*{{{*/ + + + int i,j,k; + int my_rank; + bool process_units = true; + + /*intermediary: */ + char root[50]; + int index; + int npart; + double femmodel_response; + int flag; + double *vertex_response = NULL; + double *qmu_response = NULL; + double *responses_pointer = NULL; + + /*retrieve npart: */ + parameters->FindParam(&npart,QmuNumberofpartitionsEnum); + my_rank=IssmComm::GetRank(); + + /*save the d_responses pointer: */ + responses_pointer=d_responses; + + //watch out, we have more d_numresponses than numresponsedescriptors, because the responses have been expanded if they were scaled. + //because we don't know the d_responses descriptors (the scaled ones) we can't key off them, so we will key off the responses_descriptors: */ + + for(i=0;i(vertex_response); + xDelete(qmu_response); + + } + else if (flag==IndexedEnum){ + + /*indexed response: plug index into parameters and call response module: */ + parameters->SetParam(index,IndexEnum); + + this->Responsex(&femmodel_response,root,process_units,0);//0 is the index for weights + + if(my_rank==0){ + /*plug response: */ + responses_pointer[0]=femmodel_response; + + /*increment response_pointer :*/ + responses_pointer++; + } + } + else if (flag==NodalEnum){ + _error_("nodal response functions not supported yet!"); + + /*increment response_pointer :*/ + responses_pointer++; + } + else if (flag==RegularEnum){ + + /*perfectly normal response function: */ + this->Responsex(&femmodel_response,root,process_units,0);//0 is the weight index + + if(my_rank==0){ + /*plug response: */ + responses_pointer[0]=femmodel_response; + + /*increment response_pointer :*/ + responses_pointer++; + } + } + else _error_("flag type " << flag << " not supported yet for response analysis"); + } + + /*Synthesize echo: {{{*/ + if(my_rank==0){ + _printString_(" responses: " << d_numresponses << ": "); + for(i=0;iname==MassFluxEnum){ /*to identify the mass flux that will be computed, we need the index of the profile: */ - parameters->SetParam(this->index,IndexEnum); + femmodel->parameters->SetParam(this->index,IndexEnum); } - ::Responsex(poutput_value,elements,nodes,vertices,loads,materials,parameters,this->name,false,0); + femmodel->Responsex(poutput_value,this->name,false,0); } Index: /issm/trunk-jpl/src/c/classes/objects/DependentObject.h =================================================================== --- /issm/trunk-jpl/src/c/classes/objects/DependentObject.h (revision 13698) +++ /issm/trunk-jpl/src/c/classes/objects/DependentObject.h (revision 13699) @@ -17,4 +17,5 @@ class Materials; class Parameters; +class FemModel; class DependentObject: public Object{ @@ -41,5 +42,5 @@ /*DependentObject methods: */ int NumDependents(void); - void Responsex(IssmDouble* poutput_value,Elements* elements,Nodes* nodes,Vertices* vertices,Loads* loads,Materials* materials,Parameters* parameters); + void Responsex(IssmDouble* poutput_value,FemModel* femmodel); }; Index: /issm/trunk-jpl/src/c/modules/modules.h =================================================================== --- /issm/trunk-jpl/src/c/modules/modules.h (revision 13698) +++ /issm/trunk-jpl/src/c/modules/modules.h (revision 13699) @@ -21,7 +21,5 @@ #include "./ControlInputSetGradientx/ControlInputSetGradientx.h" #include "./ControlInputScaleGradientx/ControlInputScaleGradientx.h" -#include "./CostFunctionx/CostFunctionx.h" #include "./CreateNodalConstraintsx/CreateNodalConstraintsx.h" -#include "./DakotaResponsesx/DakotaResponsesx.h" #include "./Delta18oParameterizationx/Delta18oParameterizationx.h" #include "./DragCoefficientAbsGradientx/DragCoefficientAbsGradientx.h" @@ -100,9 +98,7 @@ #include "./Reducevectorgtosx/Reducevectorgtosx.h" #include "./Reducevectorgtofx/Reducevectorgtofx.h" -#include "./RequestedOutputsx/RequestedOutputsx.h" #include "./RequestedDependentsx/RequestedDependentsx.h" #include "./ResetConstraintsx/ResetConstraintsx.h" #include "./ResetCoordinateSystemx/ResetCoordinateSystemx.h" -#include "./Responsex/Responsex.h" #include "./RheologyBbarAbsGradientx/RheologyBbarAbsGradientx.h" #include "./Scotchx/Scotchx.h" @@ -117,5 +113,4 @@ #include "./TriMeshProcessRiftsx/TriMeshProcessRiftsx.h" #include "./ThicknessAbsMisfitx/ThicknessAbsMisfitx.h" -#include "./ThicknessAbsGradientx/ThicknessAbsGradientx.h" #include "./ThicknessAlongGradientx/ThicknessAlongGradientx.h" #include "./ThicknessAcrossGradientx/ThicknessAcrossGradientx.h" Index: /issm/trunk-jpl/src/c/solutions/DakotaSpawnCore.cpp =================================================================== --- /issm/trunk-jpl/src/c/solutions/DakotaSpawnCore.cpp (revision 13698) +++ /issm/trunk-jpl/src/c/solutions/DakotaSpawnCore.cpp (revision 13699) @@ -86,5 +86,5 @@ /*compute responses: */ if(VerboseQmu()) _pprintLine_("compute dakota responses:"); - DakotaResponsesx(d_responses,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,responses_descriptors,numresponsedescriptors,d_numresponses); + femmodel->DakotaResponsesx(d_responses,responses_descriptors,numresponsedescriptors,d_numresponses); /*Free ressources:*/ Index: /issm/trunk-jpl/src/c/solutions/controltao_core.cpp =================================================================== --- /issm/trunk-jpl/src/c/solutions/controltao_core.cpp (revision 13698) +++ /issm/trunk-jpl/src/c/solutions/controltao_core.cpp (revision 13699) @@ -141,5 +141,5 @@ /*Compute objective function*/ - CostFunctionx(fcn,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters); + femmodel->CostFunctionx(fcn); /*Compute gradient*/ Index: /issm/trunk-jpl/src/c/solutions/diagnostic_core.cpp =================================================================== --- /issm/trunk-jpl/src/c/solutions/diagnostic_core.cpp (revision 13698) +++ /issm/trunk-jpl/src/c/solutions/diagnostic_core.cpp (revision 13699) @@ -102,5 +102,5 @@ InputToResultx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,PressureEnum); if(dim==3) InputToResultx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,VzEnum); - RequestedOutputsx(femmodel->results,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,requested_outputs,numoutputs); + femmodel->RequestedOutputsx(requested_outputs,numoutputs); } Index: /issm/trunk-jpl/src/c/solutions/objectivefunction.cpp =================================================================== --- /issm/trunk-jpl/src/c/solutions/objectivefunction.cpp (revision 13698) +++ /issm/trunk-jpl/src/c/solutions/objectivefunction.cpp (revision 13699) @@ -70,5 +70,5 @@ /*Compute misfit for this velocity field.*/ - CostFunctionx(&J, femmodel->elements,femmodel->nodes, femmodel->vertices,femmodel->loads, femmodel->materials, femmodel->parameters); + femmodel->CostFunctionx(&J); /*Free ressources:*/ Index: /issm/trunk-jpl/src/c/solutions/prognostic_core.cpp =================================================================== --- /issm/trunk-jpl/src/c/solutions/prognostic_core.cpp (revision 13698) +++ /issm/trunk-jpl/src/c/solutions/prognostic_core.cpp (revision 13699) @@ -51,5 +51,5 @@ if(VerboseSolution()) _pprintLine_(" saving results"); InputToResultx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,ThicknessEnum); - RequestedOutputsx(femmodel->results,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,requested_outputs,numoutputs); + femmodel->RequestedOutputsx(requested_outputs,numoutputs); } Index: /issm/trunk-jpl/src/c/solutions/steadystate_core.cpp =================================================================== --- /issm/trunk-jpl/src/c/solutions/steadystate_core.cpp (revision 13698) +++ /issm/trunk-jpl/src/c/solutions/steadystate_core.cpp (revision 13699) @@ -89,5 +89,5 @@ if(isenthalpy) InputToResultx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,EnthalpyEnum); if(!isenthalpy) InputToResultx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,BasalforcingsMeltingRateEnum); - RequestedOutputsx(femmodel->results,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,requested_outputs,numoutputs); + femmodel->RequestedOutputsx(requested_outputs,numoutputs); } Index: /issm/trunk-jpl/src/c/solutions/transient_core.cpp =================================================================== --- /issm/trunk-jpl/src/c/solutions/transient_core.cpp (revision 13698) +++ /issm/trunk-jpl/src/c/solutions/transient_core.cpp (revision 13699) @@ -140,5 +140,5 @@ InputToResultx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,SurfaceforcingsMassBalanceEnum); InputToResultx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,MaskElementonfloatingiceEnum); - RequestedOutputsx(femmodel->results,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,requested_outputs,numoutputs); + femmodel->RequestedOutputsx(requested_outputs,numoutputs); if(isdelta18o){