RheologyBAbsGradientx.cpp File Reference

#include "./RheologyBAbsGradientx.h"
#include "../../shared/shared.h"
#include "../../toolkits/toolkits.h"
#include "../../classes/Inputs2/DatasetInput2.h"

Go to the source code of this file.

Functions
void	RheologyBAbsGradientx (IssmDouble *pJ, Elements *elements, Nodes *nodes, Vertices *vertices, Loads *loads, Materials *materials, Parameters *parameters)
IssmDouble	RheologyBAbsGradient (Element *element)
void	RheologyBInitialguessMisfitx (IssmDouble *pJ, Elements *elements, Nodes *nodes, Vertices *vertices, Loads *loads, Materials *materials, Parameters *parameters)
IssmDouble	RheologyBInitialguessMisfit (Element *element)

Function Documentation

RheologyBAbsGradientx()

void RheologyBAbsGradientx	(	IssmDouble *	pJ,
		Elements *	elements,
		Nodes *	nodes,
		Vertices *	vertices,
		Loads *	loads,
		Materials *	materials,
		Parameters *	parameters
	)

Definition at line 11 of file RheologyBAbsGradientx.cpp.

                                                                                                                                                           {
 
   /*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+=RheologyBAbsGradient(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;
}

RheologyBAbsGradient()

IssmDouble RheologyBAbsGradient

(

Element *

element

)

Definition at line 32 of file RheologyBAbsGradientx.cpp.

                                                 {
 
   int         domaintype,numcomponents;
   IssmDouble  Jelem=0.;
   IssmDouble  misfit,Jdet;
   IssmDouble  dp[3],weight;
   IssmDouble* xyz_list      = NULL;
 
   /*If on water, return 0: */
   if(!element->IsIceInElement()) return 0.;
 
   /*Get problem dimension*/
   element->FindParam(&domaintype,DomainTypeEnum);
   switch(domaintype){
      case Domain2DverticalEnum:   numcomponents   = 1; break;
      case Domain3DEnum:           numcomponents   = 2; break;
      case Domain2DhorizontalEnum: numcomponents   = 2; break;
      default: _error_("not supported yet");
   }
 
   /* Get node coordinates*/
   element->GetVerticesCoordinates(&xyz_list);
 
   /*Retrieve all inputs we will be needing: */
   DatasetInput2* weights_input=element->GetDatasetInput2(InversionCostFunctionsCoefficientsEnum); _assert_(weights_input);
   Input2* rheologyb_input=element->GetInput2(MaterialsRheologyBEnum);                  _assert_(rheologyb_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,RheologyBAbsGradientEnum);
      rheologyb_input->GetInputDerivativeValue(&dp[0],xyz_list,gauss);
 
      /*Tikhonov regularization: J = 1/2 ((dp/dx)^2 + (dp/dy)^2) */ 
      Jelem+=weight*1/2*(dp[0]*dp[0] + dp[1]*dp[1])*Jdet*gauss->weight;
   }
 
   /*clean up and Return: */
   xDelete<IssmDouble>(xyz_list);
   delete gauss;
   return Jelem;
}

RheologyBInitialguessMisfitx()

void RheologyBInitialguessMisfitx	(	IssmDouble *	pJ,
		Elements *	elements,
		Nodes *	nodes,
		Vertices *	vertices,
		Loads *	loads,
		Materials *	materials,
		Parameters *	parameters
	)

Definition at line 82 of file RheologyBAbsGradientx.cpp.

                                                                                                                                                                  {
 
   /*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+=RheologyBInitialguessMisfit(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;
}

RheologyBInitialguessMisfit()

IssmDouble RheologyBInitialguessMisfit

(

Element *

element

)

Definition at line 103 of file RheologyBAbsGradientx.cpp.

                                                        {
 
   int         domaintype,numcomponents;
   IssmDouble  Jelem=0.;
   IssmDouble  misfit,Jdet;
   IssmDouble  B,B0,weight;
   IssmDouble* xyz_list      = NULL;
 
   /*If on water, return 0: */
   if(!element->IsIceInElement()) return 0.;
 
   /*Get problem dimension*/
   element->FindParam(&domaintype,DomainTypeEnum);
   switch(domaintype){
      case Domain2DverticalEnum:   numcomponents   = 1; break;
      case Domain3DEnum:           numcomponents   = 2; break;
      case Domain2DhorizontalEnum: numcomponents   = 2; break;
      default: _error_("not supported yet");
   }
 
   /* Get node coordinates*/
   element->GetVerticesCoordinates(&xyz_list);
 
   /*Retrieve all inputs we will be needing: */
   DatasetInput2* weights_input=element->GetDatasetInput2(InversionCostFunctionsCoefficientsEnum); _assert_(weights_input);
   Input2* rheologyb_input=element->GetInput2(MaterialsRheologyBbarEnum);            _assert_(rheologyb_input);
   Input2* rheologyb0_input=element->GetInput2(RheologyBInitialguessEnum);           _assert_(rheologyb0_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,RheologyBInitialguessMisfitEnum);
      rheologyb_input->GetInputValue(&B,gauss);
      rheologyb0_input->GetInputValue(&B0,gauss);
 
      /*Tikhonov regularization: J = 1/2 (B-B0)^2 */ 
      Jelem+=weight*1./2.*(B-B0)*(B-B0)*Jdet*gauss->weight;
   }
 
   /*clean up and Return: */
   xDelete<IssmDouble>(xyz_list);
   delete gauss;
   return Jelem;
}