Changeset 3180

Timestamp:

03/04/10 09:31:08 (15 years ago)

Author:

Mathieu Morlighem

Message:

Added Misfit 3

Location:

issm/trunk/src/c

Files:

• CostFunctionx/CostFunctionx.cpp (modified) (2 diffs)
• DataSet/DataSet.cpp (modified) (1 diff)
• DataSet/DataSet.h (modified) (1 diff)
• Dux/Dux.cpp (modified) (2 diffs)
• Makefile.am (modified) (2 diffs)
• Misfitx/Misfitx.cpp (modified) (2 diffs)
• SurfaceAreax (added)
• SurfaceAreax/SurfaceAreax.cpp (added)
• SurfaceAreax/SurfaceAreax.h (added)
• issm.h (modified) (1 diff)
• objects/Beam.cpp (modified) (1 diff)
• objects/Beam.h (modified) (1 diff)
• objects/Element.h (modified) (1 diff)
• objects/Input.cpp (modified) (1 diff)
• objects/Input.h (modified) (1 diff)
• objects/ParameterInputs.cpp (modified) (1 diff)
• objects/ParameterInputs.h (modified) (1 diff)
• objects/Penta.cpp (modified) (1 diff)
• objects/Penta.h (modified) (1 diff)
• objects/Sing.cpp (modified) (1 diff)
• objects/Sing.h (modified) (1 diff)
• objects/Tria.cpp (modified) (18 diffs)
• objects/Tria.h (modified) (1 diff)

issm/trunk/src/c/CostFunctionx/CostFunctionx.cpp

@@ -12 +12 @@
 #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 +28 @@
         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: */

issm/trunk/src/c/DataSet/DataSet.cpp

@@ -1530 +1530 @@
 }
 
+void  DataSet::SurfaceArea(double* pS,void* inputs,int analysis_type,int sub_analysis_type){
+
+        double S=0;;
+
+        vector<Object*>::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){
 

issm/trunk/src/c/DataSet/DataSet.h

@@ -83 +83 @@
                 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);

issm/trunk/src/c/Dux/Dux.cpp

@@ -12 +12 @@
 #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 +30 @@
         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: */

issm/trunk/src/c/Makefile.am

@@ -233 +233 @@
                                         ./ControlConstrainx/ControlConstrainx.h\
                                         ./ControlConstrainx/ControlConstrainx.cpp\
+                                        ./SurfaceAreax/SurfaceAreax.h\
+                                        ./SurfaceAreax/SurfaceAreax.cpp\
                                         ./Misfitx/Misfitx.h\
                                         ./Misfitx/Misfitx.cpp\
@@ -561 +563 @@
                                         ./ControlConstrainx/ControlConstrainx.h\
                                         ./ControlConstrainx/ControlConstrainx.cpp\
+                                        ./SurfaceAreax/SurfaceAreax.h\
+                                        ./SurfaceAreax/SurfaceAreax.cpp\
                                         ./Misfitx/Misfitx.h\
                                         ./Misfitx/Misfitx.cpp\

issm/trunk/src/c/Misfitx/Misfitx.cpp

@@ -12 +12 @@
 #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 +29 @@
         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);
 

issm/trunk/src/c/issm.h

@@ -53 +53 @@
 #include "./Orthx/Orthx.h"
 #include "./Misfitx/Misfitx.h"
+#include "./SurfaceAreax/SurfaceAreax.h"
 #include "./CostFunctionx/CostFunctionx.h"
 #include "./ControlConstrainx/ControlConstrainx.h"

issm/trunk/src/c/objects/Beam.cpp

@@ -697 +697 @@
 }
 /*}}}*/
+/*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__ 

issm/trunk/src/c/objects/Beam.h

@@ -88 +88 @@
                 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);

issm/trunk/src/c/objects/Element.h

@@ -15 +15 @@
 
         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();
                 
 };

issm/trunk/src/c/objects/Input.cpp

@@ -225 +225 @@
 }
 /*}}}*/
+/*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){

issm/trunk/src/c/objects/Input.h

@@ -53 +53 @@
                 /*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);

issm/trunk/src/c/objects/ParameterInputs.cpp

@@ -183 +183 @@
 
         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;i<dataset->Size();i++){
+                input=(Input*)dataset->GetObjectByOffset(i);
+                if (input->IsPresent(name)) return 1;
+        }
+
+        /*Input not found... return 0*/
+        return 0;
 }
 /*}}}*/

issm/trunk/src/c/objects/ParameterInputs.h

@@ -35 +35 @@
 
                 Vec  Get(char* field_name,int* dofs, int numdofs);
+                int  IsPresent(char* name);
 
                 /*declaration found in io: */

issm/trunk/src/c/objects/Penta.cpp

@@ -4209 +4209 @@
 }
 /*}}}*/
+/*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__ 

issm/trunk/src/c/objects/Penta.h

@@ -89 +89 @@
                 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);
                 

issm/trunk/src/c/objects/Sing.cpp

@@ -532 +532 @@
 }
 /*}}}*/
+/*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__ 

issm/trunk/src/c/objects/Sing.h

@@ -83 +83 @@
                 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);

issm/trunk/src/c/objects/Tria.cpp

@@ -2866 +2866 @@
         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 +2880 @@
         /* parameters: */
         double  obs_velocity_mag,velocity_mag;
-        double  absolutex,absolutey,relativex,relativey,logarithmicx,logarithmicy;
+        double  dux,duy;
 
         /*element vector : */
@@ -2900 +2896 @@
         double scaley=0;
         double scale=0;
-        double  fit=-1;
+        double S=0;
+        double fit=-1;
 
         ParameterInputs* inputs=NULL;
@@ -2918 +2915 @@
         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 +2934 @@
                 /*We are using an absolute misfit: */
                 for (i=0;i<numgrids;i++){
-                        absolutex_list[i]=obs_vx_list[i]-vx_list[i];
-                        absolutey_list[i]=obs_vy_list[i]-vy_list[i];
+                        dux_list[i]=obs_vx_list[i]-vx_list[i];
+                        duy_list[i]=obs_vy_list[i]-vy_list[i];
                 }
         }
@@ -2945 +2945 @@
                         if(obs_vx_list[i]==0)scalex=0;
                         if(obs_vy_list[i]==0)scaley=0;
-                        relativex_list[i]=scalex*(obs_vx_list[i]-vx_list[i]);
-                        relativey_list[i]=scaley*(obs_vy_list[i]-vy_list[i]);
+                        dux_list[i]=scalex*(obs_vx_list[i]-vx_list[i]);
+                        duy_list[i]=scaley*(obs_vy_list[i]-vy_list[i]);
                 }
         }
@@ -2955 +2955 @@
                         obs_velocity_mag=sqrt(pow(obs_vx_list[i],2)+pow(obs_vy_list[i],2))+numpar->epsvel; //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;i<numgrids;i++){
+                        scale=1.0/(S*sqrt(pow(vx_list[i]-obs_vx_list[i],2)+pow(vy_list[i]-obs_vx_list[i],2))+numpar->epsvel);
+                        dux_list[i]=scale*(obs_vx_list[i]-vx_list[i]);
+                        duy_list[i]=scale*(obs_vy_list[i]-vy_list[i]);
                 }
         }
@@ -2994 +3002 @@
 
                 /*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;i<numgrids;i++){
-                                due_g_gaussian[i*NDOF2+0]=absolutex*Jdet*gauss_weight*l1l2l3[i]; 
-                                due_g_gaussian[i*NDOF2+1]=absolutey*Jdet*gauss_weight*l1l2l3[i]; 
-                        }
-                }
-                else if(fit==1){
-                        /*We are using a relative misfit: */
-
-                        /*Compute relative(x/y) at gaussian point: */
-                        GetParameterValue(&relativex, &relativex_list[0],gauss_l1l2l3);
-                        GetParameterValue(&relativey, &relativey_list[0],gauss_l1l2l3);
-
-                        /*compute Du*/
-                        for (i=0;i<numgrids;i++){
-                                due_g_gaussian[i*NDOF2+0]=relativex*Jdet*gauss_weight*l1l2l3[i]; 
-                                due_g_gaussian[i*NDOF2+1]=relativey*Jdet*gauss_weight*l1l2l3[i]; 
-                        }
-                }
-                else if(fit==2){
-                        /*We are using a logarithmic misfit: */
-
-                        /*Compute logarithmic(x/y) at gaussian point: */
-                        GetParameterValue(&logarithmicx, &logarithmicx_list[0],gauss_l1l2l3);
-                        GetParameterValue(&logarithmicy, &logarithmicy_list[0],gauss_l1l2l3);
-
-                        /*compute Du*/
-                        for (i=0;i<numgrids;i++){
-                                due_g_gaussian[i*NDOF2+0]=logarithmicx*Jdet*gauss_weight*l1l2l3[i]; 
-                                due_g_gaussian[i*NDOF2+1]=logarithmicy*Jdet*gauss_weight*l1l2l3[i]; 
-                        }
-                }
-                else{
-                        /*Not supported yet! : */
-                        throw ErrorException(__FUNCT__,exprintf("%s%g","unsupported type of fit: ",fit));
+
+                /*Compute absolute(x/y) at gaussian point: */
+                GetParameterValue(&dux, &dux_list[0],gauss_l1l2l3);
+                GetParameterValue(&duy, &duy_list[0],gauss_l1l2l3);
+
+                /*compute Du*/
+                for (i=0;i<numgrids;i++){
+                        due_g_gaussian[i*NDOF2+0]=dux*Jdet*gauss_weight*l1l2l3[i]; 
+                        due_g_gaussian[i*NDOF2+1]=duy*Jdet*gauss_weight*l1l2l3[i]; 
                 }
 
@@ -3047 +3022 @@
         VecSetValues(du_g,numdof,doflist,(const double*)due_g,ADD_VALUES);
 
-cleanup_and_return: 
         xfree((void**)&first_gauss_area_coord);
         xfree((void**)&second_gauss_area_coord);
@@ -4446 +4420 @@
         double obs_vx_list[numgrids];
         double obs_vy_list[numgrids];
-        double absolute_list[numgrids];
-        double relative_list[numgrids];
-        double logarithmic_list[numgrids];
+        double misfit_list[numgrids];
 
         /* gaussian points: */
@@ -4461 +4433 @@
         /* parameters: */
         double  velocity_mag,obs_velocity_mag;
-        double  absolute,relative,logarithmic;
+        double  misfit;
 
         /* Jacobian: */
@@ -4469 +4441 @@
         double scalex=1;
         double scaley=1;
-        double  fit=-1;
+        double S=0;
+        double fit=-1;
 
         ParameterInputs* inputs=NULL;
@@ -4484 +4457 @@
         /* Recover input data: */
         if(!inputs->Recover("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 +4479 @@
                 /*We are using an absolute misfit: */
                 for (i=0;i<numgrids;i++){
-                        absolute_list[i]=0.5*(pow((vx_list[i]-obs_vx_list[i]),(double)2)+pow((vy_list[i]-obs_vy_list[i]),(double)2));
+                        misfit_list[i]=0.5*(pow((vx_list[i]-obs_vx_list[i]),(double)2)+pow((vy_list[i]-obs_vy_list[i]),(double)2));
                 }
         }
@@ -4513 +4489 @@
                         if(obs_vx_list[i]==0)scalex=0;
                         if(obs_vy_list[i]==0)scaley=0;
-                        relative_list[i]=0.5*(scalex*pow((vx_list[i]-obs_vx_list[i]),2)+scaley*pow((vy_list[i]-obs_vy_list[i]),2));
+                        misfit_list[i]=0.5*(scalex*pow((vx_list[i]-obs_vx_list[i]),2)+scaley*pow((vy_list[i]-obs_vy_list[i]),2));
                 }
         }
@@ -4521 +4497 @@
                         velocity_mag=sqrt(pow(vx_list[i],(double)2)+pow(vy_list[i],(double)2))+numpar->epsvel; //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<numgrids;i++){
+                        misfit_list[i]=sqrt(pow(vx_list[i]-obs_vx_list[i],2)+pow(vy_list[i]-obs_vx_list[i],2))/S;
                 }
         }
@@ -4543 +4525 @@
                 GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss_l1l2l3);
 
-                /*Differents misfits are allowed: */
-                if(fit==0){
-                        /*Compute absolute misfit at gaussian point: */
-                        GetParameterValue(&absolute, &absolute_list[0],gauss_l1l2l3);
-
-                        /*compute Misfit*/
-                        Jelem+=absolute*Jdet*gauss_weight;
-                }
-                else if(fit==1){
-                        /*Compute relative misfit at gaussian point: */
-                        GetParameterValue(&relative, &relative_list[0],gauss_l1l2l3);
-
-                        /*compute Misfit*/
-                        Jelem+=relative*Jdet*gauss_weight;
-                }       
-                else if(fit==2){
-                        /*Compute logarithmic misfit at gaussian point: */
-                        GetParameterValue(&logarithmic, &logarithmic_list[0],gauss_l1l2l3);
-
-                        /*compute Misfit*/
-                        Jelem+=logarithmic*Jdet*gauss_weight;
-                }
-                else throw ErrorException(__FUNCT__,exprintf("%s%i%s","fit type",fit," not supported yet!"));
-
+                /*Compute misfit at gaussian point: */
+                GetParameterValue(&misfit, &misfit_list[0],gauss_l1l2l3);
+
+                /*compute Misfit*/
+                Jelem+=misfit*Jdet*gauss_weight;
         }
 
@@ -4635 +4598 @@
         *(surface_normal+2)=normal[2]/normal_norm;
 
+}
+/*}}}*/
+/*FUNCTION SurfaceArea {{{1*/
+#undef __FUNCT__ 
+#define __FUNCT__ "Tria::SurfaceArea"
+double Tria::SurfaceArea(void* vinputs,int analysis_type,int sub_analysis_type){
+
+        int i;
+
+        /* output: */
+        double S;
+
+        /* node data: */
+        int numgrids=3;
+        double xyz_list[numgrids][3];
+        double v13[3];
+        double v23[3];
+        double normal[3];
+
+        /*If on water, return 0: */
+        if(onwater)return 0;
+
+        /* Get node coordinates and dof list: */
+        GetVerticesCoordinates(&xyz_list[0][0], nodes, numgrids);
+
+        for (i=0;i<3;i++){
+                v13[i]=xyz_list[0][i]-xyz_list[2][i];
+                v23[i]=xyz_list[1][i]-xyz_list[2][i];
+        }
+
+        normal[0]=v13[1]*v23[2]-v13[2]*v23[1];
+        normal[1]=v13[2]*v23[0]-v13[0]*v23[2];
+        normal[2]=v13[0]*v23[1]-v13[1]*v23[0];
+
+        S = 0.5 * sqrt(pow(normal[0],(double)2)+pow(normal[1],(double)2)+pow(normal[2],(double)2));
+
+        /*Return: */
+        return S;
 }
 /*}}}*/

issm/trunk/src/c/objects/Tria.h

@@ -100 +100 @@
                 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);