Changeset 17975 for issm/trunk-jpl/src/c/modules/SurfaceRelVelMisfitx/SurfaceRelVelMisfitx.cpp

Timestamp:

05/09/14 14:51:29 (11 years ago)

Author:

Mathieu Morlighem

Message:

CHG: moving misfits to modules

File:

• issm/trunk-jpl/src/c/modules/SurfaceRelVelMisfitx/SurfaceRelVelMisfitx.cpp (modified) (2 diffs)

issm/trunk-jpl/src/c/modules/SurfaceRelVelMisfitx/SurfaceRelVelMisfitx.cpp

@@ -10 +10 @@
 void SurfaceRelVelMisfitx( IssmDouble* pJ, Elements* elements,Nodes* nodes, Vertices* vertices, Loads* loads, Materials* materials,Parameters* parameters){
 
-        /*Intermediary*/
-        int i;
-        Element* element=NULL;
-
         /*output: */
-        IssmDouble J=0;
+        IssmDouble J=0.;
         IssmDouble J_sum;
 
         /*Compute Misfit: */
-        for (i=0;i<elements->Size();i++){
-                element=dynamic_cast<Element*>(elements->GetObjectByOffset(i));
-                J+=element->SurfaceRelVelMisfit();
+        for(int i=0;i<elements->Size();i++){
+                Element* element=dynamic_cast<Element*>(elements->GetObjectByOffset(i));
+                J+=SurfaceRelVelMisfit(element);
         }
 
@@ -32 +28 @@
         *pJ=J;
 }
+
+IssmDouble SurfaceRelVelMisfit(Element* element){
+
+        int        domaintype,numcomponents;
+        bool       surfaceintegral;
+        IssmDouble Jelem=0.;
+        IssmDouble misfit,Jdet,scalex,scaley;
+        IssmDouble epsvel=2.220446049250313e-16;
+        IssmDouble meanvel=3.170979198376458e-05; /*1000 m/yr*/
+        IssmDouble vx,vy,vxobs,vyobs,weight;
+        IssmDouble* xyz_list = NULL;
+
+        /*Get basal element*/
+        if(!element->IsOnSurface()) return 0.;
+
+        /*If on water, return 0: */
+        if(!element->IsIceInElement()) return 0.;
+
+        /*Get problem dimension*/
+        element->FindParam(&domaintype,DomainTypeEnum);
+        switch(domaintype){
+                case Domain2DverticalEnum:
+                        surfaceintegral = true;
+                        numcomponents   = 1;
+                        break;
+                case Domain3DEnum:
+                        surfaceintegral = true;
+                        numcomponents   = 2;
+                        break;
+                case Domain2DhorizontalEnum:
+                        surfaceintegral = false;
+                        numcomponents   = 2;
+                        break;
+                default: _error_("not supported yet");
+        }
+
+        /*Spawn surface element*/
+        Element* topelement = element->SpawnTopElement();
+
+        /* Get node coordinates*/
+        topelement->GetVerticesCoordinates(&xyz_list);
+
+        /*Retrieve all inputs we will be needing: */
+        Input* weights_input=topelement->GetInput(InversionCostFunctionsCoefficientsEnum); _assert_(weights_input);
+        Input* vx_input     =topelement->GetInput(VxEnum);                                 _assert_(vx_input);
+        Input* vxobs_input  =topelement->GetInput(InversionVxObsEnum);                     _assert_(vxobs_input);
+        Input* vy_input     = NULL;
+        Input* vyobs_input  = NULL;
+        if(numcomponents==2){
+                vy_input    =topelement->GetInput(VyEnum);              _assert_(vy_input);
+                vyobs_input =topelement->GetInput(InversionVyObsEnum);  _assert_(vyobs_input);
+        }
+
+        /* Start  looping on the number of gaussian points: */
+        Gauss* gauss=topelement->NewGauss(4);
+        for(int ig=gauss->begin();ig<gauss->end();ig++){
+
+                gauss->GaussPoint(ig);
+
+                /* Get Jacobian determinant: */
+                topelement->JacobianDeterminant(&Jdet,xyz_list,gauss);
+
+                /*Get all parameters at gaussian point*/
+                weights_input->GetInputValue(&weight,gauss,SurfaceRelVelMisfitEnum);
+                vx_input->GetInputValue(&vx,gauss);
+                vxobs_input->GetInputValue(&vxobs,gauss);
+                if(numcomponents==2){
+                        vy_input->GetInputValue(&vy,gauss);
+                        vyobs_input->GetInputValue(&vyobs,gauss);
+                }
+
+                /*Compute SurfaceRelVelMisfit:
+                 *                        
+                 *      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                      
+                 */
+
+                if(numcomponents==2){
+                        scalex=pow(meanvel/(vxobs+epsvel),2); if(vxobs==0)scalex=0;
+                        scaley=pow(meanvel/(vyobs+epsvel),2); if(vyobs==0)scaley=0;
+                        misfit=0.5*(scalex*pow((vx-vxobs),2)+scaley*pow((vy-vyobs),2));
+                }
+                else{
+                        scalex=pow(meanvel/(vxobs+epsvel),2); if(vxobs==0)scalex=0;
+                        misfit=0.5*(scalex*pow((vx-vxobs),2));
+                }
+
+                /*Add to cost function*/
+                Jelem+=misfit*weight*Jdet*gauss->weight;
+        }
+
+        /*clean up and Return: */
+        if(domaintype!=Domain2DhorizontalEnum){topelement->DeleteMaterials(); delete topelement;};
+        xDelete<IssmDouble>(xyz_list);
+        delete gauss;
+        return Jelem;
+}