ThicknessAcrossGradientx.cpp

Go to the documentation of this file.

 
#include "./ThicknessAcrossGradientx.h"
 
#include "../../shared/shared.h"
#include "../../toolkits/toolkits.h"
#include "../../classes/Inputs2/DatasetInput2.h"
 
void ThicknessAcrossGradientx( IssmDouble* pJ, Elements* elements,Nodes* nodes, Vertices* vertices, Loads* loads, Materials* materials,Parameters* parameters){
 
   /*output: */
   IssmDouble J=0.;
   IssmDouble J_sum;
 
   /*Compute Misfit: */
   for(int i=0;i<elements->Size();i++){
      Element* element=xDynamicCast<Element*>(elements->GetObjectByOffset(i));
      J+=ThicknessAcrossGradient(element);
   }
 
   /*Sum all J from all cpus of the cluster:*/
   ISSM_MPI_Reduce (&J,&J_sum,1,ISSM_MPI_DOUBLE,ISSM_MPI_SUM,0,IssmComm::GetComm() );
   ISSM_MPI_Bcast(&J_sum,1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm());
   J=J_sum;
 
   /*Assign output pointers: */
   *pJ=J;
}
 
IssmDouble ThicknessAcrossGradient(Element* element){
 
   IssmDouble  thickness,thicknessobs,weight;
   IssmDouble  Jelem=0.;
   IssmDouble  misfit,Jdet;
   IssmDouble* xyz_list = NULL;
   IssmDouble  dp[3];
   IssmDouble  vx,vy,vel;
 
   /*If on water, return 0: */
   if(!element->IsIceInElement()) return 0.;
 
   /* Get node coordinates*/
   element->GetVerticesCoordinates(&xyz_list);
 
   /*Retrieve all inputs we will be needing: */
   DatasetInput2* weights_input   =element->GetDatasetInput2(InversionCostFunctionsCoefficientsEnum); _assert_(weights_input);
   Input2* thickness_input =element->GetInput2(ThicknessEnum);                          _assert_(thickness_input);
   Input2* vx_input        =element->GetInput2(VxEnum);                                 _assert_(vx_input);
   Input2* vy_input        =element->GetInput2(VyEnum);                                 _assert_(vy_input);
 
   /* Start  looping on the number of gaussian points: */
   Gauss* gauss=element->NewGauss(2);
   for(int ig=gauss->begin();ig<gauss->end();ig++){
 
      gauss->GaussPoint(ig);
 
      /* Get Jacobian determinant: */
      element->JacobianDeterminant(&Jdet,xyz_list,gauss);
 
      /*Get all parameters at gaussian point*/
      weights_input->GetInputValue(&weight,gauss,ThicknessAcrossGradientEnum);
      thickness_input->GetInputValue(&thickness,gauss);
      thickness_input->GetInputDerivativeValue(&dp[0],xyz_list,gauss);
      vx_input->GetInputValue(&vx,gauss);
      vy_input->GetInputValue(&vy,gauss);
      vel = sqrt(vx*vx+vy*vy);
      vx  = vx/(vel+1.e-9);
      vy  = vy/(vel+1.e-9);
 
      /*J = 1/2 ( -vy*dH/dx + vx*dH/dy )^2 */
      Jelem+=weight*1/2*(-vy*dp[0] + vx*dp[1])*(-vy*dp[0] + vx*dp[1])*Jdet*gauss->weight;
   }
 
   /*clean up and Return: */
   xDelete<IssmDouble>(xyz_list);
   delete gauss;
   return Jelem;
}