BalancethicknessAnalysis Class Reference

#include <BalancethicknessAnalysis.h>

Inheritance diagram for BalancethicknessAnalysis:

Public Member Functions
void	CreateConstraints (Constraints *constraints, IoModel *iomodel)
void	CreateLoads (Loads *loads, IoModel *iomodel)
void	CreateNodes (Nodes *nodes, IoModel *iomodel, bool isamr=false)
int	DofsPerNode (int **doflist, int domaintype, int approximation)
void	UpdateElements (Elements *elements, Inputs2 *inputs2, IoModel *iomodel, int analysis_counter, int analysis_type)
void	UpdateParameters (Parameters *parameters, IoModel *iomodel, int solution_enum, int analysis_enum)
void	Core (FemModel *femmodel)
ElementVector *	CreateDVector (Element *element)
ElementMatrix *	CreateJacobianMatrix (Element *element)
ElementMatrix *	CreateKMatrix (Element *element)
ElementMatrix *	CreateKMatrixCG (Element *element)
ElementMatrix *	CreateKMatrixDG (Element *element)
ElementVector *	CreatePVector (Element *element)
ElementVector *	CreatePVectorCG (Element *element)
ElementVector *	CreatePVectorDG (Element *element)
void	GetB (IssmDouble *B, Element *element, IssmDouble *xyz_list, Gauss *gauss)
void	GetBprime (IssmDouble *B, Element *element, IssmDouble *xyz_list, Gauss *gauss)
void	GetSolutionFromInputs (Vector< IssmDouble > *solution, Element *element)
void	GradientJ (Vector< IssmDouble > *gradient, Element *element, int control_type, int control_index)
void	InputUpdateFromSolution (IssmDouble *solution, Element *element)
void	UpdateConstraints (FemModel *femmodel)
Public Member Functions inherited from Analysis
virtual	~Analysis ()

Detailed Description

Definition at line 11 of file BalancethicknessAnalysis.h.

Member Function Documentation

CreateConstraints()

void BalancethicknessAnalysis::CreateConstraints ( Constraints *  constraints, IoModel *  iomodel  )

virtual

Implements Analysis.

Definition at line 9 of file BalancethicknessAnalysis.cpp.

                                                                                         {/*{{{*/
 
   /*Fetch parameters: */
   int    stabilization;   
   iomodel->FindConstant(&stabilization,"md.balancethickness.stabilization");
 
   /*Do not add constraints in DG*/
   if(stabilization!=3){
      IoModelToConstraintsx(constraints,iomodel,"md.balancethickness.spcthickness",BalancethicknessAnalysisEnum,P1Enum);
   }
 
}/*}}}*/

CreateLoads()

void BalancethicknessAnalysis::CreateLoads ( Loads *  loads, IoModel *  iomodel  )

virtual

Implements Analysis.

Definition at line 21 of file BalancethicknessAnalysis.cpp.

                                                                        {/*{{{*/
 
   /*Intermediary*/
   int element;
   int stabilization;
 
   /*Fetch parameters: */
   iomodel->FindConstant(&stabilization,"md.balancethickness.stabilization");
 
   /*Loads only in DG*/
   if (stabilization==3){
 
      /*Get faces and elements*/
      CreateFaces(iomodel);
      iomodel->FetchData(1,"md.geometry.thickness");
 
      /*First load data:*/
      for(int i=0;i<iomodel->numberoffaces;i++){
 
         /*Get left and right elements*/
         element=iomodel->faces[4*i+2]-1; //faces are [node1 node2 elem1 elem2]
 
         /*Now, if this element is not in the partition, pass: */
         if(!iomodel->my_elements[element]) continue;
 
         /* Add load */
         loads->AddObject(new Numericalflux(i+1,i,i,iomodel));
      }
 
      /*Free data: */
      iomodel->DeleteData(1,"md.geometry.thickness");
   }
}/*}}}*/

CreateNodes()

void BalancethicknessAnalysis::CreateNodes ( Nodes *  nodes, IoModel *  iomodel, bool  isamr = false  )

virtual

Implements Analysis.

Definition at line 54 of file BalancethicknessAnalysis.cpp.

                                                                                  {/*{{{*/
 
   int  stabilization;
   iomodel->FindConstant(&stabilization,"md.balancethickness.stabilization");
 
   /*Check in 3d*/
   if(stabilization==3 && iomodel->domaintype==Domain3DEnum) _error_("DG 3d not implemented yet");
 
   /*First fetch data: */
   if(iomodel->domaintype==Domain3DEnum) iomodel->FetchData(2,"md.mesh.vertexonbase","md.mesh.vertexonsurface");
   if(stabilization!=3){
      ::CreateNodes(nodes,iomodel,BalancethicknessAnalysisEnum,P1Enum);
   }
   else{
      ::CreateNodes(nodes,iomodel,BalancethicknessAnalysisEnum,P1DGEnum);
   }
   iomodel->DeleteData(2,"md.mesh.vertexonbase","md.mesh.vertexonsurface");
}/*}}}*/

DofsPerNode()

int BalancethicknessAnalysis::DofsPerNode ( int **  doflist, int  domaintype, int  approximation  )

virtual

Implements Analysis.

Definition at line 72 of file BalancethicknessAnalysis.cpp.

                                                                                        {/*{{{*/
   return 1;
}/*}}}*/

UpdateElements()

void BalancethicknessAnalysis::UpdateElements ( Elements *  elements, Inputs2 *  inputs2, IoModel *  iomodel, int  analysis_counter, int  analysis_type  )

virtual

Implements Analysis.

Definition at line 75 of file BalancethicknessAnalysis.cpp.

                                                                                                                                        {/*{{{*/
 
   int    stabilization,finiteelement;
 
   /*Fetch data needed: */
   iomodel->FindConstant(&stabilization,"md.balancethickness.stabilization");
 
   /*Finite element type*/
   finiteelement = P1Enum;
   if(stabilization==3){
      finiteelement = P1DGEnum;
   }
 
   /*Update elements: */
   int counter=0;
   for(int i=0;i<iomodel->numberofelements;i++){
      if(iomodel->my_elements[i]){
         Element* element=(Element*)elements->GetObjectByOffset(counter);
         element->Update(inputs2,i,iomodel,analysis_counter,analysis_type,finiteelement);
         counter++;
      }
   }
 
   iomodel->FetchDataToInput(inputs2,elements,"md.geometry.thickness",ThicknessEnum);
   iomodel->FetchDataToInput(inputs2,elements,"md.geometry.surface",SurfaceEnum);
   iomodel->FetchDataToInput(inputs2,elements,"md.geometry.base",BaseEnum);
   iomodel->FetchDataToInput(inputs2,elements,"md.solidearth.sealevel",SealevelEnum,0);
   iomodel->FetchDataToInput(inputs2,elements,"md.mask.ice_levelset",MaskIceLevelsetEnum);
   iomodel->FetchDataToInput(inputs2,elements,"md.initialization.vx",VxEnum);
   iomodel->FetchDataToInput(inputs2,elements,"md.initialization.vy",VyEnum);
   iomodel->FetchDataToInput(inputs2,elements,"md.basalforcings.groundedice_melting_rate",BasalforcingsGroundediceMeltingRateEnum);
   iomodel->FetchDataToInput(inputs2,elements,"md.smb.mass_balance",SmbMassBalanceEnum);
   iomodel->FetchDataToInput(inputs2,elements,"md.balancethickness.thickening_rate",BalancethicknessThickeningRateEnum);
 
   if(iomodel->domaintype!=Domain2DhorizontalEnum){
      iomodel->FetchDataToInput(inputs2,elements,"md.mesh.vertexonbase",MeshVertexonbaseEnum);
      iomodel->FetchDataToInput(inputs2,elements,"md.mesh.vertexonsurface",MeshVertexonsurfaceEnum);
   }
}/*}}}*/

UpdateParameters()

void BalancethicknessAnalysis::UpdateParameters ( Parameters *  parameters, IoModel *  iomodel, int  solution_enum, int  analysis_enum  )

virtual

Implements Analysis.

Definition at line 114 of file BalancethicknessAnalysis.cpp.

                                                                                                                          {/*{{{*/
   parameters->AddObject(iomodel->CopyConstantObject("md.balancethickness.stabilization",BalancethicknessStabilizationEnum));
}/*}}}*/

Core()

void BalancethicknessAnalysis::Core ( FemModel *  femmodel )

virtual

Implements Analysis.

Definition at line 119 of file BalancethicknessAnalysis.cpp.

                                                               {/*{{{*/
   _error_("not implemented");
}/*}}}*/

CreateDVector()

ElementVector * BalancethicknessAnalysis::CreateDVector ( Element *  element )

virtual

Implements Analysis.

Definition at line 122 of file BalancethicknessAnalysis.cpp.

                                                                      {/*{{{*/
   /*Default, return NULL*/
   return NULL;
}/*}}}*/

CreateJacobianMatrix()

ElementMatrix * BalancethicknessAnalysis::CreateJacobianMatrix ( Element *  element )

virtual

Implements Analysis.

Definition at line 126 of file BalancethicknessAnalysis.cpp.

                                                                             {/*{{{*/
_error_("Not implemented");
}/*}}}*/

CreateKMatrix()

ElementMatrix * BalancethicknessAnalysis::CreateKMatrix ( Element *  element )

virtual

Implements Analysis.

Definition at line 129 of file BalancethicknessAnalysis.cpp.

                                                                      {/*{{{*/
 
   if(!element->IsOnBase()) return NULL;
   Element* basalelement = element->SpawnBasalElement();
 
   ElementMatrix* Ke = NULL;
   switch(element->FiniteElement()){
      case P1Enum: case P2Enum:
         Ke = CreateKMatrixCG(basalelement);
         break;
      case P1DGEnum:
         Ke = CreateKMatrixDG(basalelement);
         break;
      default:
         _error_("Element type " << EnumToStringx(element->FiniteElement()) << " not supported yet");
   }
 
   int domaintype;
   element->FindParam(&domaintype,DomainTypeEnum);
   if(domaintype!=Domain2DhorizontalEnum){basalelement->DeleteMaterials(); delete basalelement;};
   return Ke;
}/*}}}*/

CreateKMatrixCG()

ElementMatrix * BalancethicknessAnalysis::CreateKMatrixCG

(

Element *

element

)

Definition at line 151 of file BalancethicknessAnalysis.cpp.

                                                                        {/*{{{*/
 
   /*Intermediaries */
   int        stabilization;
   int        domaintype;
   IssmDouble Jdet,D_scalar,h;
   IssmDouble vel,vx,vy,dvxdx,dvydy;
   IssmDouble dvx[2],dvy[2];
   IssmDouble* xyz_list = NULL;
 
   /*Fetch number of nodes and dof for this finite element*/
   int numnodes = element->GetNumberOfNodes();
 
   /*Initialize Element vector and other vectors*/
   ElementMatrix* Ke     = element->NewElementMatrix();
   IssmDouble*    basis  = xNew<IssmDouble>(numnodes);
   IssmDouble*    dbasis = xNew<IssmDouble>(2*numnodes);
   IssmDouble     D[2][2];
 
   /*Retrieve all inputs and parameters*/
   element->GetVerticesCoordinates(&xyz_list);
   element->FindParam(&domaintype,DomainTypeEnum);
   element->FindParam(&stabilization,BalancethicknessStabilizationEnum);
   Input2* vxaverage_input=NULL;
   Input2* vyaverage_input=NULL;
   if(domaintype==Domain2DhorizontalEnum){
      vxaverage_input=element->GetInput2(VxEnum); _assert_(vxaverage_input);
      vyaverage_input=element->GetInput2(VyEnum); _assert_(vyaverage_input);
   }
   else{
      vxaverage_input=element->GetInput2(VxAverageEnum); _assert_(vxaverage_input);
      vyaverage_input=element->GetInput2(VyAverageEnum); _assert_(vyaverage_input);
   }
   h = element->CharacteristicLength();
 
   /* 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);
 
      element->JacobianDeterminant(&Jdet,xyz_list,gauss);
      element->NodalFunctions(basis,gauss);
      element->NodalFunctionsDerivatives(dbasis,xyz_list,gauss);
 
      vxaverage_input->GetInputValue(&vx,gauss);
      vyaverage_input->GetInputValue(&vy,gauss);
      vxaverage_input->GetInputDerivativeValue(&dvx[0],xyz_list,gauss);
      vyaverage_input->GetInputDerivativeValue(&dvy[0],xyz_list,gauss);
      dvxdx=dvx[0];
      dvydy=dvy[1];
      D_scalar=gauss->weight*Jdet;
 
      for(int i=0;i<numnodes;i++){
         for(int j=0;j<numnodes;j++){
            /*\phi_i \phi_j \nabla\cdot v*/
            Ke->values[i*numnodes+j] += D_scalar*basis[i]*basis[j]*(dvxdx+dvydy);
            /*\phi_i v\cdot\nabla\phi_j*/
            Ke->values[i*numnodes+j] += D_scalar*basis[i]*(vx*dbasis[0*numnodes+j] + vy*dbasis[1*numnodes+j]);
         }
      }
 
      if(stabilization==1){
         /*Streamline upwinding*/
         vel=sqrt(vx*vx+vy*vy);
         D[0][0]=h/(2*vel)*vx*vx;
         D[1][0]=h/(2*vel)*vy*vx;
         D[0][1]=h/(2*vel)*vx*vy;
         D[1][1]=h/(2*vel)*vy*vy;
      }
      else if(stabilization==2){
         /*SSA*/
         vxaverage_input->GetInputAverage(&vx);
         vyaverage_input->GetInputAverage(&vy);
         D[0][0]=h/2.0*fabs(vx);
         D[0][1]=0.;
         D[1][0]=0.;
         D[1][1]=h/2.0*fabs(vy);
      }
      if(stabilization==1 || stabilization==2){
         D[0][0]=D_scalar*D[0][0];
         D[1][0]=D_scalar*D[1][0];
         D[0][1]=D_scalar*D[0][1];
         D[1][1]=D_scalar*D[1][1];
 
         for(int i=0;i<numnodes;i++){
            for(int j=0;j<numnodes;j++){
               Ke->values[i*numnodes+j] += (
                        dbasis[0*numnodes+i] *(D[0][0]*dbasis[0*numnodes+j] + D[0][1]*dbasis[1*numnodes+j]) +
                        dbasis[1*numnodes+i] *(D[1][0]*dbasis[0*numnodes+j] + D[1][1]*dbasis[1*numnodes+j]) 
                        );
            }
         }
      }
   }
 
   /*Clean up and return*/
   xDelete<IssmDouble>(xyz_list);
   xDelete<IssmDouble>(basis);
   xDelete<IssmDouble>(dbasis);
   delete gauss;
   return Ke;
}/*}}}*/

CreateKMatrixDG()

ElementMatrix * BalancethicknessAnalysis::CreateKMatrixDG

(

Element *

element

)

Definition at line 253 of file BalancethicknessAnalysis.cpp.

                                                                        {/*{{{*/
 
   /*Intermediaries */
   int        domaintype;
   IssmDouble Jdet,D_scalar,vx,vy,dvxdx,dvydy,vel;
   IssmDouble dvx[2],dvy[2];
   IssmDouble* xyz_list = NULL;
 
   /*Fetch number of nodes and dof for this finite element*/
   int numnodes = element->GetNumberOfNodes();
 
   /*Initialize Element vector and other vectors*/
   ElementMatrix* Ke     = element->NewElementMatrix();
   IssmDouble*    B      = xNew<IssmDouble>(2*numnodes);
   IssmDouble*    Bprime = xNew<IssmDouble>(2*numnodes);
   IssmDouble     D[2][2];
 
   /*Retrieve all inputs and parameters*/
   element->GetVerticesCoordinates(&xyz_list);
   element->FindParam(&domaintype,DomainTypeEnum);
   Input2* vxaverage_input=NULL;
   Input2* vyaverage_input=NULL;
   if(domaintype==Domain2DhorizontalEnum){
      vxaverage_input=element->GetInput2(VxEnum); _assert_(vxaverage_input);
      vyaverage_input=element->GetInput2(VyEnum); _assert_(vyaverage_input);
   }
   else{
      vxaverage_input=element->GetInput2(VxAverageEnum); _assert_(vxaverage_input);
      vyaverage_input=element->GetInput2(VyAverageEnum); _assert_(vyaverage_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);
 
      element->JacobianDeterminant(&Jdet,xyz_list,gauss);
      vxaverage_input->GetInputValue(&vx,gauss);
      vyaverage_input->GetInputValue(&vy,gauss);
      vxaverage_input->GetInputDerivativeValue(&dvx[0],xyz_list,gauss);
      vyaverage_input->GetInputDerivativeValue(&dvy[0],xyz_list,gauss);
      D_scalar=gauss->weight*Jdet;
 
      /*WARNING: B and Bprime are inverted compared to CG*/
      GetB(Bprime,element,xyz_list,gauss);
      GetBprime(B,element,xyz_list,gauss);
 
      D_scalar = - gauss->weight*Jdet;
      D[0][0]  = D_scalar*vx;
      D[0][1]  = 0.;
      D[1][0]  = 0.;
      D[1][1]  = D_scalar*vy;
      TripleMultiply(B,2,numnodes,1,
               &D[0][0],2,2,0,
               Bprime,2,numnodes,0,
               &Ke->values[0],1);
 
   }
 
   /*Clean up and return*/
   xDelete<IssmDouble>(xyz_list);
   xDelete<IssmDouble>(B);
   xDelete<IssmDouble>(Bprime);
   delete gauss;
   return Ke;
}/*}}}*/

CreatePVector()

ElementVector * BalancethicknessAnalysis::CreatePVector ( Element *  element )

virtual

Implements Analysis.

Definition at line 319 of file BalancethicknessAnalysis.cpp.

                                                                      {/*{{{*/
 
   if(!element->IsOnBase()) return NULL;
   Element* basalelement = element->SpawnBasalElement();
 
   ElementVector* pe = NULL;
   switch(element->FiniteElement()){
      case P1Enum: case P2Enum:
         pe = CreatePVectorCG(basalelement);
         break;
      case P1DGEnum:
         pe = CreatePVectorDG(basalelement);
         break;
      default:
         _error_("Element type " << EnumToStringx(element->FiniteElement()) << " not supported yet");
   }
 
   int domaintype;
   element->FindParam(&domaintype,DomainTypeEnum);
   if(domaintype!=Domain2DhorizontalEnum){basalelement->DeleteMaterials(); delete basalelement;};
   return pe;
}/*}}}*/

CreatePVectorCG()

ElementVector * BalancethicknessAnalysis::CreatePVectorCG

(

Element *

element

)

Definition at line 341 of file BalancethicknessAnalysis.cpp.

                                                                        {/*{{{*/
 
   /*Intermediaries */
   IssmDouble  dhdt,mb,ms,Jdet;
   IssmDouble* xyz_list = NULL;
 
   /*Fetch number of nodes and dof for this finite element*/
   int numnodes = element->GetNumberOfNodes();
 
   /*Initialize Element vector and other vectors*/
   ElementVector* pe    = element->NewElementVector();
   IssmDouble*    basis = xNew<IssmDouble>(numnodes);
 
   /*Retrieve all inputs and parameters*/
   element->GetVerticesCoordinates(&xyz_list);
   Input2* mb_input   = element->GetInput2(BasalforcingsGroundediceMeltingRateEnum);       _assert_(mb_input);
   Input2* ms_input   = element->GetInput2(SmbMassBalanceEnum);     _assert_(ms_input);
   Input2* dhdt_input = element->GetInput2(BalancethicknessThickeningRateEnum); _assert_(dhdt_input);
 
   /*Initialize mb_correction to 0, do not forget!:*/
   /* 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);
 
      element->JacobianDeterminant(&Jdet,xyz_list,gauss);
      element->NodalFunctions(basis,gauss);
 
      ms_input->GetInputValue(&ms,gauss);
      mb_input->GetInputValue(&mb,gauss);
      dhdt_input->GetInputValue(&dhdt,gauss);
 
      for(int i=0;i<numnodes;i++) pe->values[i]+=Jdet*gauss->weight*(ms-mb-dhdt)*basis[i];
   }
 
   /*Clean up and return*/
   xDelete<IssmDouble>(xyz_list);
   xDelete<IssmDouble>(basis);
   delete gauss;
   return pe;
}/*}}}*/

CreatePVectorDG()

ElementVector * BalancethicknessAnalysis::CreatePVectorDG

(

Element *

element

)

Definition at line 382 of file BalancethicknessAnalysis.cpp.

                                                                        {/*{{{*/
 
   /*Intermediaries */
   IssmDouble  dhdt,mb,ms,Jdet;
   IssmDouble* xyz_list = NULL;
 
   /*Fetch number of nodes and dof for this finite element*/
   int numnodes = element->GetNumberOfNodes();
 
   /*Initialize Element vector and other vectors*/
   ElementVector* pe    = element->NewElementVector();
   IssmDouble*    basis = xNew<IssmDouble>(numnodes);
 
   /*Retrieve all inputs and parameters*/
   element->GetVerticesCoordinates(&xyz_list);
   Input2* mb_input   = element->GetInput2(BasalforcingsGroundediceMeltingRateEnum);       _assert_(mb_input);
   Input2* ms_input   = element->GetInput2(SmbMassBalanceEnum);     _assert_(ms_input);
   Input2* dhdt_input = element->GetInput2(BalancethicknessThickeningRateEnum); _assert_(dhdt_input);
 
   /*Initialize mb_correction to 0, do not forget!:*/
   /* 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);
 
      element->JacobianDeterminant(&Jdet,xyz_list,gauss);
      element->NodalFunctions(basis,gauss);
 
      ms_input->GetInputValue(&ms,gauss);
      mb_input->GetInputValue(&mb,gauss);
      dhdt_input->GetInputValue(&dhdt,gauss);
 
      for(int i=0;i<numnodes;i++) pe->values[i]+=Jdet*gauss->weight*(ms-mb-dhdt)*basis[i];
   }
 
   /*Clean up and return*/
   xDelete<IssmDouble>(xyz_list);
   xDelete<IssmDouble>(basis);
   delete gauss;
   return pe;
}/*}}}*/

GetB()

void BalancethicknessAnalysis::GetB	(	IssmDouble *	B,
		Element *	element,
		IssmDouble *	xyz_list,
		Gauss *	gauss
	)

Definition at line 423 of file BalancethicknessAnalysis.cpp.

                                                                                                             {/*{{{*/
   /*Compute B  matrix. B=[B1 B2 B3] where Bi is of size 3*2. 
    * For node i, Bi can be expressed in the actual coordinate system
    * by: 
    *       Bi=[ N ]
    *          [ N ]
    * where N is the finiteelement function for node i.
    *
    * We assume B_prog has been allocated already, of size: 2x(1*numnodes)
    */
 
   /*Fetch number of nodes for this finite element*/
   int numnodes = element->GetNumberOfNodes();
 
   /*Get nodal functions*/
   IssmDouble* basis=xNew<IssmDouble>(numnodes);
   element->NodalFunctions(basis,gauss);
 
   /*Build B: */
   for(int i=0;i<numnodes;i++){
      B[numnodes*0+i] = basis[i];
      B[numnodes*1+i] = basis[i];
   }
 
   /*Clean-up*/
   xDelete<IssmDouble>(basis);
}/*}}}*/

GetBprime()

void BalancethicknessAnalysis::GetBprime	(	IssmDouble *	B,
		Element *	element,
		IssmDouble *	xyz_list,
		Gauss *	gauss
	)

Definition at line 450 of file BalancethicknessAnalysis.cpp.

                                                                                                                       {/*{{{*/
   /*Compute B'  matrix. B'=[B1' B2' B3'] where Bi' is of size 3*2. 
    * For node i, Bi' can be expressed in the actual coordinate system
    * by: 
    *       Bi_prime=[ dN/dx ]
    *                [ dN/dy ]
    * where N is the finiteelement function for node i.
    *
    * We assume B' has been allocated already, of size: 3x(2*numnodes)
    */
 
   /*Fetch number of nodes for this finite element*/
   int numnodes = element->GetNumberOfNodes();
 
   /*Get nodal functions derivatives*/
   IssmDouble* dbasis=xNew<IssmDouble>(2*numnodes);
   element->NodalFunctionsDerivatives(dbasis,xyz_list,gauss);
 
   /*Build B': */
   for(int i=0;i<numnodes;i++){
      Bprime[numnodes*0+i] = dbasis[0*numnodes+i];
      Bprime[numnodes*1+i] = dbasis[1*numnodes+i];
   }
 
   /*Clean-up*/
   xDelete<IssmDouble>(dbasis);
 
}/*}}}*/

GetSolutionFromInputs()

void BalancethicknessAnalysis::GetSolutionFromInputs ( Vector< IssmDouble > *  solution, Element *  element  )

virtual

Implements Analysis.

Definition at line 478 of file BalancethicknessAnalysis.cpp.

                                                                                                           {/*{{{*/
      _error_("not implemented yet");
}/*}}}*/

GradientJ()

void BalancethicknessAnalysis::GradientJ ( Vector< IssmDouble > *  gradient, Element *  element, int  control_type, int  control_index  )

virtual

Implements Analysis.

Definition at line 481 of file BalancethicknessAnalysis.cpp.

                                                                                                                                  {/*{{{*/
   /* WARNING: this gradient is valid for Soft balance thickness only */
 
   /*If on water, grad = 0: */
   if(!element->IsIceInElement()) return;
 
   /*Intermediaries*/
   IssmDouble Jdet,weight;
   IssmDouble thickness,thicknessobs,dH[3],dp[3];
   IssmDouble  vx,vy,vel,dvx[2],dvy[2],dhdt,basal_melting,surface_mass_balance;
   IssmDouble *xyz_list= NULL;
 
   /*Get list of cost functions*/
   int *responses = NULL;
   int  num_responses,resp,solution;
   element->FindParam(&num_responses,InversionNumCostFunctionsEnum);
   element->FindParam(&responses,NULL,InversionCostFunctionsEnum);
   element->FindParam(&solution,SolutionTypeEnum);
   if(solution!=BalancethicknessSoftSolutionEnum) _error_("not implemented yet");
   if(control_type!=ThicknessEnum)                _error_("Control "<<EnumToStringx(control_type)<<" not supported");
 
   /*Fetch number of vertices for this finite element*/
   int numvertices = element->GetNumberOfVertices();
 
   /*Initialize some vectors*/
   IssmDouble* basis         = xNew<IssmDouble>(numvertices);
   IssmDouble* dbasis        = xNew<IssmDouble>(2*numvertices);
   IssmDouble* ge            = xNewZeroInit<IssmDouble>(numvertices);
   int*        vertexpidlist = xNew<int>(numvertices);
 
   /*Retrieve all inputs we will be needing: */
   element->GetVerticesCoordinates(&xyz_list);
   element->GradientIndexing(&vertexpidlist[0],control_index);
   DatasetInput2* weights_input       = element->GetDatasetInput2(InversionCostFunctionsCoefficientsEnum);  _assert_(weights_input);
   Input2* thickness_input            = element->GetInput2(ThicknessEnum);                           _assert_(thickness_input);
   Input2* thicknessobs_input         = element->GetInput2(InversionThicknessObsEnum);               _assert_(thicknessobs_input);
   Input2* vx_input                   = element->GetInput2(VxEnum);                                  _assert_(vx_input);
   Input2* vy_input                   = element->GetInput2(VyEnum);                                  _assert_(vy_input);
   Input2* surface_mass_balance_input = element->GetInput2(SmbMassBalanceEnum);          _assert_(surface_mass_balance_input);
   Input2* basal_melting_input        = element->GetInput2(BasalforcingsGroundediceMeltingRateEnum); _assert_(basal_melting_input);
   Input2* dhdt_input                 = element->GetInput2(BalancethicknessThickeningRateEnum);      _assert_(dhdt_input);
 
   /* Start  looping on the number of gaussian points: */
   Gauss* gauss=element->NewGauss(4);
   for(int ig=gauss->begin();ig<gauss->end();ig++){
      gauss->GaussPoint(ig);
 
      thickness_input->GetInputValue(&thickness, gauss);
      thickness_input->GetInputDerivativeValue(&dH[0],xyz_list,gauss);
      thicknessobs_input->GetInputValue(&thicknessobs, gauss);
 
      element->JacobianDeterminant(&Jdet,xyz_list,gauss);
      element->NodalFunctionsP1(basis,gauss);
      element->NodalFunctionsP1Derivatives(dbasis,xyz_list,gauss);
 
      /*Deal with first part (partial derivative a J with respect to k)*/
      for(resp=0;resp<num_responses;resp++){
 
         weights_input->GetInputValue(&weight,gauss,responses[resp]);
 
         switch(responses[resp]){
            case ThicknessAbsMisfitEnum:
               for(int i=0;i<numvertices;i++) ge[i]+= (thicknessobs-thickness)*weight*Jdet*gauss->weight*basis[i];
               break;
            case ThicknessAbsGradientEnum:
               for(int i=0;i<numvertices;i++) ge[i]+= - weight*dH[0]*dbasis[0*numvertices+i]*Jdet*gauss->weight;
               for(int i=0;i<numvertices;i++) ge[i]+= - weight*dH[1]*dbasis[1*numvertices+i]*Jdet*gauss->weight;
               break;
            case ThicknessAlongGradientEnum:
               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);
               for(int i=0;i<numvertices;i++) ge[i]+= - weight*(dH[0]*vx+dH[1]*vy)*(dbasis[0*numvertices+i]*vx+dbasis[1*numvertices+i]*vy)*Jdet*gauss->weight;
               break;
            case ThicknessAcrossGradientEnum:
               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);
               for(int i=0;i<numvertices;i++) ge[i]+= - weight*(dH[0]*(-vy)+dH[1]*vx)*(dbasis[0*numvertices+i]*(-vy)+dbasis[1*numvertices+i]*vx)*Jdet*gauss->weight;
               break;
            case BalancethicknessMisfitEnum:
               surface_mass_balance_input->GetInputValue(&surface_mass_balance,gauss);
               basal_melting_input->GetInputValue(&basal_melting,gauss);
               dhdt_input->GetInputValue(&dhdt,gauss);
               vx_input->GetInputValue(&vx,gauss);
               vx_input->GetInputDerivativeValue(&dvx[0],xyz_list,gauss);
               vy_input->GetInputValue(&vy,gauss);
               vy_input->GetInputDerivativeValue(&dvy[0],xyz_list,gauss);
               for(int i=0;i<numvertices;i++){
                  ge[i]+= - weight*Jdet*gauss->weight*(
                     (vx*dH[0]+vy*dH[1] + thickness*(dvx[0]+dvy[1]))*(vx*dbasis[0*numvertices+i]+ vy*dbasis[1*numvertices+i] + basis[i]*(dvx[0]+dvy[1]))
                     -(surface_mass_balance-basal_melting-dhdt)*(vx*dbasis[0*numvertices+i]+ vy*dbasis[1*numvertices+i] + basis[i]*(dvx[0]+dvy[1]))
                     );
               }
               break;
            default:
               _error_("response " << EnumToStringx(responses[resp]) << " not supported yet");
         }
      }
   }
   gradient->SetValues(numvertices,vertexpidlist,ge,ADD_VAL);
 
   /*Clean up and return*/
   xDelete<IssmDouble>(xyz_list);
   xDelete<IssmDouble>(basis);
   xDelete<IssmDouble>(ge);
   xDelete<int>(vertexpidlist);
   xDelete<int>(responses);
   delete gauss;
 
}/*}}}*/

InputUpdateFromSolution()

void BalancethicknessAnalysis::InputUpdateFromSolution ( IssmDouble *  solution, Element *  element  )

virtual

Implements Analysis.

Definition at line 596 of file BalancethicknessAnalysis.cpp.

                                                                                                     {/*{{{*/
 
   int domaintype;
   element->FindParam(&domaintype,DomainTypeEnum);
   switch(domaintype){
      case Domain2DhorizontalEnum:
         element->InputUpdateFromSolutionOneDof(solution,ThicknessEnum);
         break;
      case Domain3DEnum:
         element->InputUpdateFromSolutionOneDofCollapsed(solution,ThicknessEnum);
         break;
      default: _error_("mesh "<<EnumToStringx(domaintype)<<" not supported yet");
   }
}/*}}}*/

UpdateConstraints()

void BalancethicknessAnalysis::UpdateConstraints ( FemModel *  femmodel )

virtual

Implements Analysis.

Definition at line 610 of file BalancethicknessAnalysis.cpp.

                                                                            {/*{{{*/
   /*Default, do nothing*/
   return;
}/*}}}*/

---

The documentation for this class was generated from the following files:

• src/c/analyses/BalancethicknessAnalysis.h
• src/c/analyses/BalancethicknessAnalysis.cpp