Index: /issm/trunk/src/c/CostFunctionx/CostFunctionx.cpp =================================================================== --- /issm/trunk/src/c/CostFunctionx/CostFunctionx.cpp (revision 3179) +++ /issm/trunk/src/c/CostFunctionx/CostFunctionx.cpp (revision 3180) @@ -12,8 +12,13 @@ #include "../toolkits/toolkits.h" #include "../EnumDefinitions/EnumDefinitions.h" +#include "../SurfaceAreax/SurfaceAreax.h" void CostFunctionx( double* pJ, DataSet* elements,DataSet* nodes, DataSet* loads, DataSet* materials,DataSet* parameters, ParameterInputs* inputs,int analysis_type,int sub_analysis_type){ - + + /*Intermediary*/ + double fit; + double S; + /*output: */ double J; @@ -23,4 +28,12 @@ elements->Configure(elements,loads, nodes, materials,parameters); parameters->Configure(elements,loads, nodes, materials,parameters); + + /*If fit=3, compute Surface Area*/ + if(!inputs->Recover("fit",&fit)) throw ErrorException(__FUNCT__," missing fit input parameter"); + if (fit==3 && !inputs->IsPresent("surfacearea")){ + + SurfaceAreax(&S,elements,nodes,loads,materials,parameters,inputs,analysis_type,sub_analysis_type); + inputs->Add("surfacearea",S); + } /*Compute gradients: */ Index: /issm/trunk/src/c/DataSet/DataSet.cpp =================================================================== --- /issm/trunk/src/c/DataSet/DataSet.cpp (revision 3179) +++ /issm/trunk/src/c/DataSet/DataSet.cpp (revision 3180) @@ -1530,4 +1530,26 @@ } +void DataSet::SurfaceArea(double* pS,void* inputs,int analysis_type,int sub_analysis_type){ + + double S=0;; + + vector::iterator object; + Element* element=NULL; + + for ( object=objects.begin() ; object < objects.end(); object++ ){ + + if(EnumIsElement((*object)->Enum())){ + + element=(Element*)(*object); + S+=element->SurfaceArea(inputs,analysis_type,sub_analysis_type); + + } + } + + /*Assign output pointers:*/ + *pS=S; + +} + 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 3179) +++ /issm/trunk/src/c/DataSet/DataSet.h (revision 3180) @@ -83,4 +83,5 @@ 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 SurfaceArea(double* pS, 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/Dux/Dux.cpp =================================================================== --- /issm/trunk/src/c/Dux/Dux.cpp (revision 3179) +++ /issm/trunk/src/c/Dux/Dux.cpp (revision 3180) @@ -12,11 +12,15 @@ #include "../toolkits/toolkits.h" #include "../EnumDefinitions/EnumDefinitions.h" +#include "../SurfaceAreax/SurfaceAreax.h" void Dux( Vec* pdu_g, DataSet* elements,DataSet* nodes, DataSet* loads, DataSet* materials, DataSet* parameters, ParameterInputs* inputs,int analysis_type,int sub_analysis_type){ + /*Intermediary*/ int i; int gsize; int found; + double fit=-1; + double S; /*output: */ @@ -26,4 +30,12 @@ elements->Configure(elements,loads, nodes, materials,parameters); parameters->Configure(elements,loads, nodes, materials,parameters); + + /*If fit=3, compute Surface Area*/ + if(!inputs->Recover("fit",&fit)) throw ErrorException(__FUNCT__," missing fit input parameter"); + if (fit==3 && !inputs->IsPresent("surfacearea")){ + + SurfaceAreax(&S,elements,nodes,loads,materials,parameters,inputs,analysis_type,sub_analysis_type); + inputs->Add("surfacearea",S); + } /*Get size of matrix: */ Index: /issm/trunk/src/c/Makefile.am =================================================================== --- /issm/trunk/src/c/Makefile.am (revision 3179) +++ /issm/trunk/src/c/Makefile.am (revision 3180) @@ -233,4 +233,6 @@ ./ControlConstrainx/ControlConstrainx.h\ ./ControlConstrainx/ControlConstrainx.cpp\ + ./SurfaceAreax/SurfaceAreax.h\ + ./SurfaceAreax/SurfaceAreax.cpp\ ./Misfitx/Misfitx.h\ ./Misfitx/Misfitx.cpp\ @@ -561,4 +563,6 @@ ./ControlConstrainx/ControlConstrainx.h\ ./ControlConstrainx/ControlConstrainx.cpp\ + ./SurfaceAreax/SurfaceAreax.h\ + ./SurfaceAreax/SurfaceAreax.cpp\ ./Misfitx/Misfitx.h\ ./Misfitx/Misfitx.cpp\ Index: /issm/trunk/src/c/Misfitx/Misfitx.cpp =================================================================== --- /issm/trunk/src/c/Misfitx/Misfitx.cpp (revision 3179) +++ /issm/trunk/src/c/Misfitx/Misfitx.cpp (revision 3180) @@ -12,8 +12,13 @@ #include "../toolkits/toolkits.h" #include "../EnumDefinitions/EnumDefinitions.h" +#include "../SurfaceAreax/SurfaceAreax.h" void Misfitx( double* pJ, DataSet* elements,DataSet* nodes, DataSet* loads, DataSet* materials,DataSet* parameters, ParameterInputs* inputs,int analysis_type,int sub_analysis_type){ + /*Intermediary*/ + int fit; + double S; + /*output: */ double J; @@ -24,5 +29,13 @@ parameters->Configure(elements,loads, nodes, materials,parameters); - /*Compute gradients: */ + /*If fit=3, compute Surface Area*/ + if(!inputs->Recover("fit",&fit)) throw ErrorException(__FUNCT__," missing fit input parameter"); + if (fit==3 && !inputs->IsPresent("surfacearea")){ + + SurfaceAreax(&S,elements,nodes,loads,materials,parameters,inputs,analysis_type,sub_analysis_type); + inputs->Add("surfacearea",S); + } + + /*Compute Misfit: */ elements->Misfit(&J,inputs,analysis_type,sub_analysis_type); Index: /issm/trunk/src/c/SurfaceAreax/SurfaceAreax.cpp =================================================================== --- /issm/trunk/src/c/SurfaceAreax/SurfaceAreax.cpp (revision 3180) +++ /issm/trunk/src/c/SurfaceAreax/SurfaceAreax.cpp (revision 3180) @@ -0,0 +1,36 @@ +/*!\file SurfaceAreax + * \brief: compute Surface area + */ + +#include "./SurfaceAreax.h" + +#undef __FUNCT__ +#define __FUNCT__ "SurfaceAreax" + +#include "../shared/shared.h" +#include "../include/macros.h" +#include "../toolkits/toolkits.h" +#include "../EnumDefinitions/EnumDefinitions.h" + +void SurfaceAreax( double* pS, DataSet* elements,DataSet* nodes, DataSet* loads, DataSet* materials,DataSet* parameters, + ParameterInputs* inputs,int analysis_type,int sub_analysis_type){ + + /*output: */ + double S; + double S_sum; + + /*First, get elements and loads configured: */ + elements->Configure(elements,loads, nodes, materials,parameters); + parameters->Configure(elements,loads, nodes, materials,parameters); + + /*Compute gradients: */ + elements->SurfaceArea(&S,inputs,analysis_type,sub_analysis_type); + + /*Sum all J from all cpus of the cluster:*/ + MPI_Reduce (&S,&S_sum,1,MPI_DOUBLE,MPI_SUM,0,MPI_COMM_WORLD ); + MPI_Bcast(&S_sum,1,MPI_DOUBLE,0,MPI_COMM_WORLD); + S=S_sum; + + /*Assign output pointers: */ + *pS=S; +} Index: /issm/trunk/src/c/SurfaceAreax/SurfaceAreax.h =================================================================== --- /issm/trunk/src/c/SurfaceAreax/SurfaceAreax.h (revision 3180) +++ /issm/trunk/src/c/SurfaceAreax/SurfaceAreax.h (revision 3180) @@ -0,0 +1,15 @@ +/*!\file: SurfaceAreax.h + * \brief header file for ... + */ + +#ifndef _SURFACEAREAX_H +#define _SURFACEAREAX_H + +#include "../DataSet/DataSet.h" + +/* local prototypes: */ +void SurfaceAreax( double* pS, DataSet* elements,DataSet* nodes, DataSet* loads, DataSet* materials, DataSet* parameters, + ParameterInputs* inputs,int analysis_type,int sub_analysis_type); + +#endif /* _SURFACEAREAX_H */ + Index: /issm/trunk/src/c/issm.h =================================================================== --- /issm/trunk/src/c/issm.h (revision 3179) +++ /issm/trunk/src/c/issm.h (revision 3180) @@ -53,4 +53,5 @@ #include "./Orthx/Orthx.h" #include "./Misfitx/Misfitx.h" +#include "./SurfaceAreax/SurfaceAreax.h" #include "./CostFunctionx/CostFunctionx.h" #include "./ControlConstrainx/ControlConstrainx.h" Index: /issm/trunk/src/c/objects/Beam.cpp =================================================================== --- /issm/trunk/src/c/objects/Beam.cpp (revision 3179) +++ /issm/trunk/src/c/objects/Beam.cpp (revision 3180) @@ -697,4 +697,11 @@ } /*}}}*/ +/*FUNCTION Beam SurfaceArea{{{1*/ +#undef __FUNCT__ +#define __FUNCT__ "Beam::SurfaceArea" +double Beam::SurfaceArea(void*,int,int){ + throw ErrorException(__FUNCT__," not supported yet!"); +} +/*}}}*/ /*FUNCTION Beam UpdateFromInputs{{{1*/ #undef __FUNCT__ Index: /issm/trunk/src/c/objects/Beam.h =================================================================== --- /issm/trunk/src/c/objects/Beam.h (revision 3179) +++ /issm/trunk/src/c/objects/Beam.h (revision 3180) @@ -88,4 +88,5 @@ void GradjB(_p_Vec*, void*, int,int ); double Misfit(void*,int,int); + double SurfaceArea(void*,int,int); double CostFunction(void*,int,int); void GetNodalFunctions(double* l1l2, double gauss_coord); Index: /issm/trunk/src/c/objects/Element.h =================================================================== --- /issm/trunk/src/c/objects/Element.h (revision 3179) +++ /issm/trunk/src/c/objects/Element.h (revision 3180) @@ -15,33 +15,34 @@ public: - virtual ~Element(){}; - virtual void Echo()=0; - virtual void DeepEcho()=0; - virtual int GetId()=0; - virtual int MyRank()=0; - virtual void Marshall(char** pmarshalled_dataset)=0; - virtual int MarshallSize()=0; - virtual char* GetName()=0; - virtual void Demarshall(char** pmarshalled_dataset)=0; - virtual void Configure(void* loads,void* nodes,void* materials,void* parameters)=0; - virtual void CreateKMatrix(Mat Kgg,void* inputs,int analysis_type,int sub_analysis_type)=0; - virtual void CreatePVector(Vec pg, void* inputs, int analysis_type,int sub_analysis_type)=0; - virtual void UpdateFromInputs(void* inputs)=0; - virtual void GetNodes(void** nodes)=0; - virtual void* GetMatPar()=0; - virtual int GetShelf()=0; - virtual int GetOnBed()=0; - virtual void GetThicknessList(double* thickness_list)=0; - virtual void GetBedList(double* bed_list)=0; - virtual void Du(Vec du_g,void* inputs,int analysis_type,int sub_analysis_type)=0; - virtual void Gradj(Vec grad_g,void* inputs,int analysis_type,int sub_analysis_type,char* control_type)=0; - 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 ~Element(){}; + virtual void Echo()=0; + virtual void DeepEcho()=0; + virtual int GetId()=0; + virtual int MyRank()=0; + virtual void Marshall(char** pmarshalled_dataset)=0; + virtual int MarshallSize()=0; + virtual char* GetName()=0; + virtual void Demarshall(char** pmarshalled_dataset)=0; + virtual void Configure(void* loads,void* nodes,void* materials,void* parameters)=0; + virtual void CreateKMatrix(Mat Kgg,void* inputs,int analysis_type,int sub_analysis_type)=0; + virtual void CreatePVector(Vec pg, void* inputs, int analysis_type,int sub_analysis_type)=0; + virtual void UpdateFromInputs(void* inputs)=0; + virtual void GetNodes(void** nodes)=0; + virtual void* GetMatPar()=0; + virtual int GetShelf()=0; + virtual int GetOnBed()=0; + virtual void GetThicknessList(double* thickness_list)=0; + virtual void GetBedList(double* bed_list)=0; + virtual void Du(Vec du_g,void* inputs,int analysis_type,int sub_analysis_type)=0; + virtual void Gradj(Vec grad_g,void* inputs,int analysis_type,int sub_analysis_type,char* control_type)=0; + 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)=0; virtual double CostFunction(void* inputs,int analysis_type,int sub_analysis_type)=0; - virtual void ComputePressure(Vec p_g)=0; + virtual double SurfaceArea(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; - int Enum(); + int Enum(); }; Index: /issm/trunk/src/c/objects/Input.cpp =================================================================== --- /issm/trunk/src/c/objects/Input.cpp (revision 3179) +++ /issm/trunk/src/c/objects/Input.cpp (revision 3180) @@ -225,4 +225,10 @@ } /*}}}*/ +/*FUNCTION Input::IsPresent {{{1*/ +int Input::IsPresent(char* input_name){ + if (strcmp(input_name,name)==0) return 1; + else return 0; +} +/*}}}*/ /*FUNCTION Input::IsSameName {{{1*/ int Input::IsSameName(char* input_name){ Index: /issm/trunk/src/c/objects/Input.h =================================================================== --- /issm/trunk/src/c/objects/Input.h (revision 3179) +++ /issm/trunk/src/c/objects/Input.h (revision 3180) @@ -53,4 +53,5 @@ /*Input specific routines: */ int IsSameName(char* input_name); + int IsPresent(char* input_name); void Recover(double* values, int ndof, int* dofs,int numnodes,void** nodes); void Recover(int* pinteger); Index: /issm/trunk/src/c/objects/ParameterInputs.cpp =================================================================== --- /issm/trunk/src/c/objects/ParameterInputs.cpp (revision 3179) +++ /issm/trunk/src/c/objects/ParameterInputs.cpp (revision 3180) @@ -183,4 +183,21 @@ return ug; +} +/*}}}*/ +/*FUNCTION ParameterInputs::IsPresent(char* name,char* string) {{{1*/ +int ParameterInputs::IsPresent(char* name){ + + /*Intermediary*/ + int i; + Input* input=NULL; + + /*Go through dataset, and figure out if the input is present*/ + for(i=0;iSize();i++){ + input=(Input*)dataset->GetObjectByOffset(i); + if (input->IsPresent(name)) return 1; + } + + /*Input not found... return 0*/ + return 0; } /*}}}*/ Index: /issm/trunk/src/c/objects/ParameterInputs.h =================================================================== --- /issm/trunk/src/c/objects/ParameterInputs.h (revision 3179) +++ /issm/trunk/src/c/objects/ParameterInputs.h (revision 3180) @@ -35,4 +35,5 @@ Vec Get(char* field_name,int* dofs, int numdofs); + int IsPresent(char* name); /*declaration found in io: */ Index: /issm/trunk/src/c/objects/Penta.cpp =================================================================== --- /issm/trunk/src/c/objects/Penta.cpp (revision 3179) +++ /issm/trunk/src/c/objects/Penta.cpp (revision 3180) @@ -4209,4 +4209,39 @@ } /*}}}*/ +/*FUNCTION SurfaceArea {{{1*/ +#undef __FUNCT__ +#define __FUNCT__ "Penta::SurfaceArea" +double Penta::SurfaceArea(void* inputs,int analysis_type,int sub_analysis_type){ + + double S; + 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 SurfaceArea*/ + tria=(Tria*)SpawnTria(0,1,2); //grids 0, 1 and 2 make the new tria (lower face). + S=tria->SurfaceArea(inputs,analysis_type,sub_analysis_type); + delete tria; + return S; + } + else{ + + tria=(Tria*)SpawnTria(3,4,5); //grids 3, 4 and 5 make the new tria (upper face). + S=tria->SurfaceArea(inputs,analysis_type,sub_analysis_type); + delete tria; + return S; + } +} +/*}}}*/ /*FUNCTION SurfaceNormal {{{1*/ #undef __FUNCT__ Index: /issm/trunk/src/c/objects/Penta.h =================================================================== --- /issm/trunk/src/c/objects/Penta.h (revision 3179) +++ /issm/trunk/src/c/objects/Penta.h (revision 3180) @@ -89,4 +89,5 @@ 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); + double SurfaceArea(void* inputs,int analysis_type,int sub_analysis_type); double CostFunction(void* inputs,int analysis_type,int sub_analysis_type); Index: /issm/trunk/src/c/objects/Sing.cpp =================================================================== --- /issm/trunk/src/c/objects/Sing.cpp (revision 3179) +++ /issm/trunk/src/c/objects/Sing.cpp (revision 3180) @@ -532,4 +532,11 @@ } /*}}}*/ +/*FUNCTION Sing::SurfaceArea {{{1*/ +#undef __FUNCT__ +#define __FUNCT__ "Sing::SurfaceArea" +double Sing::SurfaceArea(void*, int,int){ + throw ErrorException(__FUNCT__," not supported yet!"); +} +/*}}}*/ /*FUNCTION Sing::UpdateFromInputs {{{1*/ #undef __FUNCT__ Index: /issm/trunk/src/c/objects/Sing.h =================================================================== --- /issm/trunk/src/c/objects/Sing.h (revision 3179) +++ /issm/trunk/src/c/objects/Sing.h (revision 3180) @@ -83,4 +83,5 @@ void GradjB(_p_Vec*, void*, int,int); double Misfit(void*,int,int); + double SurfaceArea(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 3179) +++ /issm/trunk/src/c/objects/Tria.cpp (revision 3180) @@ -2866,10 +2866,6 @@ double obs_vx_list[numgrids]; double obs_vy_list[numgrids]; - double absolutex_list[numgrids]; - double absolutey_list[numgrids]; - double relativex_list[numgrids]; - double relativey_list[numgrids]; - double logarithmicx_list[numgrids]; - double logarithmicy_list[numgrids]; + double dux_list[numgrids]; + double duy_list[numgrids]; /* gaussian points: */ @@ -2884,5 +2880,5 @@ /* parameters: */ double obs_velocity_mag,velocity_mag; - double absolutex,absolutey,relativex,relativey,logarithmicx,logarithmicy; + double dux,duy; /*element vector : */ @@ -2900,5 +2896,6 @@ double scaley=0; double scale=0; - double fit=-1; + double S=0; + double fit=-1; ParameterInputs* inputs=NULL; @@ -2918,4 +2915,7 @@ if(!inputs->Recover("velocity_obs",&obs_vxvy_list[0][0],2,dofs2,numgrids,(void**)nodes)){ throw ErrorException(__FUNCT__,"missing velocity_obs input parameter"); + } + if(fit==3 && !inputs->Recover("surfacearea",&S)){ + throw ErrorException(__FUNCT__,"missing surface area input parameter"); } if(!inputs->Recover("velocity",&vxvy_list[0][0],2,dofs2,numgrids,(void**)nodes)){ @@ -2934,6 +2934,6 @@ /*We are using an absolute misfit: */ for (i=0;iepsvel; //epsvel to avoid observed velocity being nil. scale=-8*pow(numpar->meanvel,2)/pow(velocity_mag,2)*log(velocity_mag/obs_velocity_mag); - logarithmicx_list[i]=scale*vx_list[i]; - logarithmicy_list[i]=scale*vy_list[i]; + dux_list[i]=scale*vx_list[i]; + duy_list[i]=scale*vy_list[i]; + } + } + else if(fit==3){ + /*We are using an spacially average absolute misfit: */ + for (i=0;iepsvel); + dux_list[i]=scale*(obs_vx_list[i]-vx_list[i]); + duy_list[i]=scale*(obs_vy_list[i]-vy_list[i]); } } @@ -2994,46 +3002,13 @@ /*Build due_g_gaussian vector: we have three cases here, according to which type of misfit we are using. */ - if(fit==0){ - /*We are using an absolute misfit: */ - - /*Compute absolute(x/y) at gaussian point: */ - GetParameterValue(&absolutex, &absolutex_list[0],gauss_l1l2l3); - GetParameterValue(&absolutey, &absolutey_list[0],gauss_l1l2l3); - - /*compute Du*/ - for (i=0;iRecover("fit",&fit)) throw ErrorException(__FUNCT__," missing fit input parameter"); + if(fit==3 && !inputs->Recover("surfacearea",&S)){ + throw ErrorException(__FUNCT__,"missing surface area input parameter"); + } if(!inputs->Recover("velocity_obs",&obs_vxvy_list[0][0],2,dofs2,numgrids,(void**)nodes)){ throw ErrorException(__FUNCT__,"missing velocity_obs input parameter"); @@ -4503,5 +4479,5 @@ /*We are using an absolute misfit: */ for (i=0;iepsvel; //epsvel to avoid velocity being nil. obs_velocity_mag=sqrt(pow(obs_vx_list[i],(double)2)+pow(obs_vy_list[i],(double)2))+numpar->epsvel; //epsvel to avoid observed velocity being nil. - logarithmic_list[i]=4*pow(numpar->meanvel,(double)2)*pow(log(velocity_mag/obs_velocity_mag),(double)2); + misfit_list[i]=4*pow(numpar->meanvel,(double)2)*pow(log(velocity_mag/obs_velocity_mag),(double)2); + } + } + else if(fit==3){ + /*We are using an spacially average absolute misfit: */ + for (i=0;i