Index: /issm/trunk/src/c/CostFunctionx/CostFunctionx.cpp =================================================================== --- /issm/trunk/src/c/CostFunctionx/CostFunctionx.cpp (revision 3169) +++ /issm/trunk/src/c/CostFunctionx/CostFunctionx.cpp (revision 3169) @@ -0,0 +1,36 @@ +/*!\file CostFunctionx + * \brief: compute misfit between observations and model + */ + +#include "./CostFunctionx.h" + +#undef __FUNCT__ +#define __FUNCT__ "CostFunctionx" + +#include "../shared/shared.h" +#include "../include/macros.h" +#include "../toolkits/toolkits.h" +#include "../EnumDefinitions/EnumDefinitions.h" + +void CostFunctionx( double* pJ, DataSet* elements,DataSet* nodes, DataSet* loads, DataSet* materials,DataSet* parameters, + ParameterInputs* inputs,int analysis_type,int sub_analysis_type){ + + /*output: */ + double J; + double J_sum; + + /*First, get elements and loads configured: */ + elements->Configure(elements,loads, nodes, materials,parameters); + parameters->Configure(elements,loads, nodes, materials,parameters); + + /*Compute gradients: */ + elements->CostFunction(&J,inputs,analysis_type,sub_analysis_type); + + /*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; + + /*Assign output pointers: */ + *pJ=J; +} Index: /issm/trunk/src/c/CostFunctionx/CostFunctionx.h =================================================================== --- /issm/trunk/src/c/CostFunctionx/CostFunctionx.h (revision 3169) +++ /issm/trunk/src/c/CostFunctionx/CostFunctionx.h (revision 3169) @@ -0,0 +1,15 @@ +/*!\file: CostFunctionx.h + * \brief header file for inverse methods misfit computation + */ + +#ifndef _COSTFUNCTIONX_H +#define _COSTFUNCTIONX_H + +#include "../DataSet/DataSet.h" + +/* local prototypes: */ +void CostFunctionx( double* pJ, DataSet* elements,DataSet* nodes, DataSet* loads, DataSet* materials, DataSet* parameters, + ParameterInputs* inputs,int analysis_type,int sub_analysis_type); + +#endif /* _MISFITX_H */ + Index: /issm/trunk/src/c/DataSet/DataSet.cpp =================================================================== --- /issm/trunk/src/c/DataSet/DataSet.cpp (revision 3168) +++ /issm/trunk/src/c/DataSet/DataSet.cpp (revision 3169) @@ -1486,5 +1486,5 @@ } -void DataSet::Misfit(double* pJ,void* inputs,int analysis_type,int sub_analysis_type,int real){ +void DataSet::Misfit(double* pJ,void* inputs,int analysis_type,int sub_analysis_type){ double J=0;; @@ -1498,5 +1498,5 @@ element=(Element*)(*object); - J+=element->Misfit(inputs,analysis_type,sub_analysis_type,real); + J+=element->Misfit(inputs,analysis_type,sub_analysis_type); } @@ -1508,4 +1508,26 @@ } +void DataSet::CostFunction(double* pJ,void* inputs,int analysis_type,int sub_analysis_type){ + + double J=0;; + + vector::iterator object; + Element* element=NULL; + + for ( object=objects.begin() ; object < objects.end(); object++ ){ + + if(EnumIsElement((*object)->Enum())){ + + element=(Element*)(*object); + J+=element->CostFunction(inputs,analysis_type,sub_analysis_type); + + } + } + + /*Assign output pointers:*/ + *pJ=J; + +} + void DataSet::FieldExtrude(Vec field,double* field_serial,char* field_name, int collapse){ Index: /issm/trunk/src/c/DataSet/DataSet.h =================================================================== --- /issm/trunk/src/c/DataSet/DataSet.h (revision 3168) +++ /issm/trunk/src/c/DataSet/DataSet.h (revision 3169) @@ -81,5 +81,6 @@ void Du(Vec du_g,void* inputs,int analysis_type,int sub_analysis_type); void Gradj(Vec grad_g,void* inputs,int analysis_type,int sub_analysis_type,char* control_type); - void Misfit(double* pJ, void* inputs,int analysis_type,int sub_analysis_type,int real); + void Misfit(double* pJ, void* inputs,int analysis_type,int sub_analysis_type); + void CostFunction(double* pJ, void* inputs,int analysis_type,int sub_analysis_type); void FieldDepthAverageAtBase(Vec field,double* field_serial,char* fieldname); int DeleteObject(Object* object); Index: /issm/trunk/src/c/Makefile.am =================================================================== --- /issm/trunk/src/c/Makefile.am (revision 3168) +++ /issm/trunk/src/c/Makefile.am (revision 3169) @@ -235,4 +235,6 @@ ./Misfitx/Misfitx.h\ ./Misfitx/Misfitx.cpp\ + ./CostFunctionx/CostFunctionx.h\ + ./CostFunctionx/CostFunctionx.cpp\ ./Orthx/Orthx.h\ ./Orthx/Orthx.cpp\ @@ -561,4 +563,6 @@ ./Misfitx/Misfitx.h\ ./Misfitx/Misfitx.cpp\ + ./CostFunctionx/CostFunctionx.h\ + ./CostFunctionx/CostFunctionx.cpp\ ./Orthx/Orthx.h\ ./Orthx/Orthx.cpp\ Index: /issm/trunk/src/c/Misfitx/Misfitx.cpp =================================================================== --- /issm/trunk/src/c/Misfitx/Misfitx.cpp (revision 3168) +++ /issm/trunk/src/c/Misfitx/Misfitx.cpp (revision 3169) @@ -14,5 +14,5 @@ void Misfitx( double* pJ, DataSet* elements,DataSet* nodes, DataSet* loads, DataSet* materials,DataSet* parameters, - ParameterInputs* inputs,int analysis_type,int sub_analysis_type,int real){ + ParameterInputs* inputs,int analysis_type,int sub_analysis_type){ /*output: */ @@ -25,5 +25,5 @@ /*Compute gradients: */ - elements->Misfit(&J,inputs,analysis_type,sub_analysis_type,real); + elements->Misfit(&J,inputs,analysis_type,sub_analysis_type); /*Sum all J from all cpus of the cluster:*/ Index: /issm/trunk/src/c/Misfitx/Misfitx.h =================================================================== --- /issm/trunk/src/c/Misfitx/Misfitx.h (revision 3168) +++ /issm/trunk/src/c/Misfitx/Misfitx.h (revision 3169) @@ -10,5 +10,5 @@ /* local prototypes: */ void Misfitx( double* pJ, DataSet* elements,DataSet* nodes, DataSet* loads, DataSet* materials, DataSet* parameters, - ParameterInputs* inputs,int analysis_type,int sub_analysis_type,int real); + ParameterInputs* inputs,int analysis_type,int sub_analysis_type); #endif /* _MISFITX_H */ Index: /issm/trunk/src/c/Qmux/DakotaResponses.cpp =================================================================== --- /issm/trunk/src/c/Qmux/DakotaResponses.cpp (revision 3168) +++ /issm/trunk/src/c/Qmux/DakotaResponses.cpp (revision 3169) @@ -264,5 +264,5 @@ /*Compute misfit: */ - Misfitx( &J, femmodel->elements,femmodel->nodes, femmodel->loads, femmodel->materials, femmodel->parameters,inputs,analysis_type,sub_analysis_type,1); + Misfitx( &J, femmodel->elements,femmodel->nodes, femmodel->loads, femmodel->materials, femmodel->parameters,inputs,analysis_type,sub_analysis_type); Index: /issm/trunk/src/c/issm.h =================================================================== --- /issm/trunk/src/c/issm.h (revision 3168) +++ /issm/trunk/src/c/issm.h (revision 3169) @@ -53,4 +53,5 @@ #include "./Orthx/Orthx.h" #include "./Misfitx/Misfitx.h" +#include "./CostFunctionx/CostFunctionx.h" #include "./ControlConstrainx/ControlConstrainx.h" #include "./FieldDepthAveragex/FieldDepthAveragex.h" Index: /issm/trunk/src/c/objects/Beam.cpp =================================================================== --- /issm/trunk/src/c/objects/Beam.cpp (revision 3168) +++ /issm/trunk/src/c/objects/Beam.cpp (revision 3169) @@ -233,4 +233,11 @@ } /*}}}*/ +/*FUNCTION Beam CostFunction{{{1*/ +#undef __FUNCT__ +#define __FUNCT__ "Beam::CostFunction" +double Beam::CostFunction(void*,int,int){ + throw ErrorException(__FUNCT__," not supported yet!"); +} +/*}}}*/ /*FUNCTION Beam CreateKMatrix{{{1*/ #undef __FUNCT__ @@ -666,5 +673,5 @@ #undef __FUNCT__ #define __FUNCT__ "Beam::Misfit" -double Beam::Misfit(void*,int,int,int){ +double Beam::Misfit(void*,int,int){ throw ErrorException(__FUNCT__," not supported yet!"); } Index: /issm/trunk/src/c/objects/Beam.h =================================================================== --- /issm/trunk/src/c/objects/Beam.h (revision 3168) +++ /issm/trunk/src/c/objects/Beam.h (revision 3169) @@ -87,5 +87,6 @@ void GradjDrag(_p_Vec*, void*, int,int ); void GradjB(_p_Vec*, void*, int,int ); - double Misfit(void*,int,int,int); + double Misfit(void*,int,int); + double CostFunction(void*,int,int); void GetNodalFunctions(double* l1l2, double gauss_coord); void GetParameterValue(double* pvalue, double* value_list,double gauss_coord); Index: /issm/trunk/src/c/objects/Element.h =================================================================== --- /issm/trunk/src/c/objects/Element.h (revision 3168) +++ /issm/trunk/src/c/objects/Element.h (revision 3169) @@ -38,5 +38,6 @@ virtual void GradjDrag(Vec grad_g,void* inputs,int analysis_type,int sub_analysis_type)=0; virtual void GradjB(Vec grad_g,void* inputs,int analysis_type,int sub_analysis_type)=0; - virtual double Misfit(void* inputs,int analysis_type,int sub_analysis_type,int real)=0; + virtual double Misfit(void* inputs,int analysis_type,int sub_analysis_type)=0; + virtual double CostFunction(void* inputs,int analysis_type,int sub_analysis_type)=0; virtual void ComputePressure(Vec p_g)=0; virtual double MassFlux(double* segment,double* ug)=0; Index: /issm/trunk/src/c/objects/Penta.cpp =================================================================== --- /issm/trunk/src/c/objects/Penta.cpp (revision 3168) +++ /issm/trunk/src/c/objects/Penta.cpp (revision 3169) @@ -282,4 +282,39 @@ Object* Penta::copy() { return new Penta(*this); +} +/*}}}*/ +/*FUNCTION CostFunction {{{1*/ +#undef __FUNCT__ +#define __FUNCT__ "Penta::CostFunction" +double Penta::CostFunction(void* inputs,int analysis_type,int sub_analysis_type){ + + double J; + Tria* tria=NULL; + + /*If on water, return 0: */ + if(onwater)return 0; + + /*Bail out if this element if: + * -> Non collapsed and not on the surface + * -> collapsed (2d model) and not on bed) */ + if ((!collapse && !onsurface) || (collapse && !onbed)){ + return 0; + } + else if (collapse){ + + /*This element should be collapsed into a tria element at its base. Create this tria element, + * and compute CostFunction*/ + tria=(Tria*)SpawnTria(0,1,2); //grids 0, 1 and 2 make the new tria (lower face). + J=tria->CostFunction(inputs,analysis_type,sub_analysis_type); + delete tria; + return J; + } + else{ + + tria=(Tria*)SpawnTria(3,4,5); //grids 3, 4 and 5 make the new tria (upper face). + J=tria->CostFunction(inputs,analysis_type,sub_analysis_type); + delete tria; + return J; + } } /*}}}*/ @@ -3915,5 +3950,5 @@ #undef __FUNCT__ #define __FUNCT__ "Penta::Misfit" -double Penta::Misfit(void* inputs,int analysis_type,int sub_analysis_type,int real){ +double Penta::Misfit(void* inputs,int analysis_type,int sub_analysis_type){ double J; @@ -3934,5 +3969,5 @@ * and compute Misfit*/ tria=(Tria*)SpawnTria(0,1,2); //grids 0, 1 and 2 make the new tria (lower face). - J=tria->Misfit(inputs,analysis_type,sub_analysis_type,real); + J=tria->Misfit(inputs,analysis_type,sub_analysis_type); delete tria; return J; @@ -3941,5 +3976,5 @@ tria=(Tria*)SpawnTria(3,4,5); //grids 3, 4 and 5 make the new tria (upper face). - J=tria->Misfit(inputs,analysis_type,sub_analysis_type,real); + J=tria->Misfit(inputs,analysis_type,sub_analysis_type); delete tria; return J; Index: /issm/trunk/src/c/objects/Penta.h =================================================================== --- /issm/trunk/src/c/objects/Penta.h (revision 3168) +++ /issm/trunk/src/c/objects/Penta.h (revision 3169) @@ -88,5 +88,6 @@ void GradjDrag(Vec grad_g,void* inputs,int analysis_type,int sub_analysis_type); void GradjB(Vec grad_g,void* inputs,int analysis_type,int sub_analysis_type); - double Misfit(void* inputs,int analysis_type,int sub_analysis_type,int real); + double Misfit(void* inputs,int analysis_type,int sub_analysis_type); + double CostFunction(void* inputs,int analysis_type,int sub_analysis_type); void GetThicknessList(double* thickness_list); Index: /issm/trunk/src/c/objects/Sing.cpp =================================================================== --- /issm/trunk/src/c/objects/Sing.cpp (revision 3168) +++ /issm/trunk/src/c/objects/Sing.cpp (revision 3169) @@ -216,4 +216,11 @@ } /*}}}*/ +/*FUNCTION Sing::CostFunction {{{1*/ +#undef __FUNCT__ +#define __FUNCT__ "Sing::CostFunction" +double Sing::CostFunction(void*, int,int){ + throw ErrorException(__FUNCT__," not supported yet!"); +} +/*}}}*/ /*FUNCTION Sing::CreateKMatrix {{{1*/ #undef __FUNCT__ @@ -515,5 +522,5 @@ #undef __FUNCT__ #define __FUNCT__ "Sing::Misfit" -double Sing::Misfit(void*, int,int,int){ +double Sing::Misfit(void*, int,int){ throw ErrorException(__FUNCT__," not supported yet!"); } Index: /issm/trunk/src/c/objects/Sing.h =================================================================== --- /issm/trunk/src/c/objects/Sing.h (revision 3168) +++ /issm/trunk/src/c/objects/Sing.h (revision 3169) @@ -82,5 +82,6 @@ void GradjDrag(_p_Vec*, void*, int,int); void GradjB(_p_Vec*, void*, int,int); - double Misfit(void*,int,int,int); + double Misfit(void*,int,int); + double CostFunction(void*,int,int); double MassFlux(double* segment,double* ug); Index: /issm/trunk/src/c/objects/Tria.cpp =================================================================== --- /issm/trunk/src/c/objects/Tria.cpp (revision 3168) +++ /issm/trunk/src/c/objects/Tria.cpp (revision 3169) @@ -268,4 +268,110 @@ } /*}}}*/ +/*FUNCTION CostFunction {{{1*/ +#undef __FUNCT__ +#define __FUNCT__ "Tria::CostFunction" +double Tria::CostFunction(void* vinputs,int analysis_type,int sub_analysis_type){ + + int i; + + /* output: */ + double Jelem; + + /* node data: */ + const int numgrids=3; + const int numdof=2*numgrids; + const int NDOF2=2; + int dofs1[1]={0}; + int dofs2[2]={0,1}; + double xyz_list[numgrids][3]; + + /* grid data: */ + double B[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 k_gauss; + double B_gauss; + + /* parameters: */ + double dk[NDOF2]; + double dB[NDOF2]; + + /* Jacobian: */ + double Jdet; + + ParameterInputs* inputs=NULL; + + /*If on water, return 0: */ + if(onwater)return 0; + + /*recover pointers: */ + inputs=(ParameterInputs*)vinputs; + + /* Get node coordinates and dof list: */ + GetVerticesCoordinates(&xyz_list[0][0], nodes, numgrids); + + /*First, get Misfit*/ + Jelem=Misfit(inputs,analysis_type,sub_analysis_type); + + /* Get gaussian points and weights (make this a statically initialized list of points? fstd): */ + GaussTria( &num_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &gauss_weights, 2); + +#ifdef _ISSM_DEBUG_ + for (i=0;icontrol_type,"drag")==0){ + if (!shelf){ + + GetParameterDerivativeValue(&dk[0], &k[0],&xyz_list[0][0], gauss_l1l2l3); + Jelem+=numpar->cm_noisedmp*1/2*(pow(dk[0],2)+pow(dk[1],2))*Jdet*gauss_weight; + + } + } + else if (strcmp(numpar->control_type,"B")==0){ + + if(!inputs->Recover("B",&B[0],1,dofs1,numgrids,(void**)nodes)){ + throw ErrorException(__FUNCT__,"parameter B not found in input"); + } + GetParameterDerivativeValue(&dB[0], &B[0],&xyz_list[0][0], gauss_l1l2l3); + Jelem+=numpar->cm_noisedmp*1/2*(pow(dB[0],2)+pow(dB[1],2))*Jdet*gauss_weight; + + } + else{ + throw ErrorException(__FUNCT__,exprintf("%s%s","unsupported control type: ",numpar->control_type)); + } + + } + + xfree((void**)&first_gauss_area_coord); + xfree((void**)&second_gauss_area_coord); + xfree((void**)&third_gauss_area_coord); + xfree((void**)&gauss_weights); + + /*Return: */ + return Jelem; +} +/*}}}*/ /*FUNCTION CreateKMatrix {{{1*/ #undef __FUNCT__ @@ -4318,5 +4424,5 @@ #undef __FUNCT__ #define __FUNCT__ "Tria::Misfit" -double Tria::Misfit(void* vinputs,int analysis_type,int sub_analysis_type,int real){ +double Tria::Misfit(void* vinputs,int analysis_type,int sub_analysis_type){ int i; @@ -4343,5 +4449,4 @@ double relative_list[numgrids]; double logarithmic_list[numgrids]; - double B[numgrids]; /* gaussian points: */ @@ -4353,12 +4458,8 @@ double gauss_weight; double gauss_l1l2l3[3]; - double k_gauss; - double B_gauss; /* parameters: */ double velocity_mag,obs_velocity_mag; double absolute,relative,logarithmic; - double dk[NDOF2]; - double dB[NDOF2]; /* Jacobian: */ @@ -4431,10 +4532,4 @@ GaussTria( &num_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &gauss_weights, 2); -#ifdef _ISSM_DEBUG_ - for (i=0;icontrol_type,"drag")==0){ - if (!shelf){ - - //noise dampening - GetParameterDerivativeValue(&dk[0], &k[0],&xyz_list[0][0], gauss_l1l2l3); - Jelem+=numpar->cm_noisedmp*1/2*(pow(dk[0],2)+pow(dk[1],2))*Jdet*gauss_weight; - - } - } - else if (strcmp(numpar->control_type,"B")==0){ - if(!inputs->Recover("B",&B[0],1,dofs1,numgrids,(void**)nodes)){ - throw ErrorException(__FUNCT__,"parameter B not found in input"); - } - //noise dampening - GetParameterDerivativeValue(&dB[0], &B[0],&xyz_list[0][0], gauss_l1l2l3); - Jelem+=numpar->cm_noisedmp*1/2*(pow(dB[0],2)+pow(dB[1],2))*Jdet*gauss_weight; - - //min dampening - GetParameterValue(&B_gauss, &B[0],gauss_l1l2l3); - if(B_gausscm_mindmp_value){ - Jelem+=numpar->cm_mindmp_slope*B_gauss*Jdet*gauss_weight; - } - - //max dampening - if(B_gauss>numpar->cm_maxdmp_value){ - Jelem+=numpar->cm_maxdmp_slope*B_gauss*Jdet*gauss_weight; - } - } - else{ - throw ErrorException(__FUNCT__,exprintf("%s%s","unsupported control type: ",numpar->control_type)); - } - } /*Differents misfits are allowed: */ @@ -4509,8 +4566,7 @@ } else throw ErrorException(__FUNCT__,exprintf("%s%i%s","fit type",fit," not supported yet!")); - - - } -cleanup_and_return: + + } + xfree((void**)&first_gauss_area_coord); xfree((void**)&second_gauss_area_coord); Index: /issm/trunk/src/c/objects/Tria.h =================================================================== --- /issm/trunk/src/c/objects/Tria.h (revision 3168) +++ /issm/trunk/src/c/objects/Tria.h (revision 3169) @@ -99,5 +99,6 @@ void SurfaceNormal(double* surface_normal, double xyz_list[3][3]); void GradjB(Vec grad_g,void* inputs,int analysis_type,int sub_analysis_type); - double Misfit(void* inputs,int analysis_type,int sub_analysis_type,int real); + double Misfit(void* inputs,int analysis_type,int sub_analysis_type); + double CostFunction(void* inputs,int analysis_type,int sub_analysis_type); void CreatePVectorDiagnosticHoriz(Vec pg,void* inputs,int analysis_type,int sub_analysis_type); Index: /issm/trunk/src/c/parallel/objectivefunctionC.cpp =================================================================== --- /issm/trunk/src/c/parallel/objectivefunctionC.cpp (revision 3168) +++ /issm/trunk/src/c/parallel/objectivefunctionC.cpp (revision 3169) @@ -97,5 +97,5 @@ /*Compute misfit for this velocity field.*/ inputs->Add("fit",fit[n]); - Misfitx( &J, femmodel->elements,femmodel->nodes, femmodel->loads, femmodel->materials, femmodel->parameters,inputs,analysis_type,sub_analysis_type,0); + CostFunctionx( &J, femmodel->elements,femmodel->nodes, femmodel->loads, femmodel->materials, femmodel->parameters,inputs,analysis_type,sub_analysis_type); /*Free ressources:*/