StressbalanceVerticalAnalysis Class Reference

#include <StressbalanceVerticalAnalysis.h>

Inheritance diagram for StressbalanceVerticalAnalysis:

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 *	CreateKMatrixBase (Element *element)
ElementMatrix *	CreateKMatrixSurface (Element *element)
ElementMatrix *	CreateKMatrixVolume (Element *element)
ElementVector *	CreatePVector (Element *element)
ElementVector *	CreatePVectorBase (Element *element)
ElementVector *	CreatePVectorSurface (Element *element)
ElementVector *	CreatePVectorVolume (Element *element)
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 StressbalanceVerticalAnalysis.h.

Member Function Documentation

CreateConstraints()

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

virtual

Implements Analysis.

Definition at line 9 of file StressbalanceVerticalAnalysis.cpp.

                                                                                              {/*{{{*/
 
   /*Intermediary*/
   bool        isSIA,isSSA,isL1L2,isHO,isFS,iscoupling;
   int         Mz,Nz;
   IssmDouble *spcvz = NULL;
 
   /*return if not 3d mesh*/
   if(iomodel->domaintype!=Domain3DEnum) return;
 
   /*Fetch parameters: */
   iomodel->FindConstant(&isSIA,"md.flowequation.isSIA");
   iomodel->FindConstant(&isSSA,"md.flowequation.isSSA");
   iomodel->FindConstant(&isL1L2,"md.flowequation.isL1L2");
   iomodel->FindConstant(&isHO,"md.flowequation.isHO");
   iomodel->FindConstant(&isFS,"md.flowequation.isFS");
 
   /*Do we have coupling*/
   if((isSIA?1.:0.) + (isSSA?1.:0.) + (isL1L2?1.:0.) + (isHO?1.:0.) + (isFS?1.:0.) >1.)
    iscoupling = true;
   else
    iscoupling = false;
 
   /*If no coupling, call Regular IoModelToConstraintsx, else, use P1 elements only*/
   if(!iscoupling){
      IoModelToConstraintsx(constraints,iomodel,"md.stressbalance.spcvz",StressbalanceVerticalAnalysisEnum,P1Enum,0);
   }
   else{
      /*Fetch data: */
      iomodel->FetchData(1,"md.flowequation.borderFS");
      /*Fetch Spc*/
      iomodel->FetchData(&spcvz,&Mz,&Nz,"md.stressbalance.spcvz");
      if(Nz>1) _error_("not supported yet (needs to be coded)");
 
      /*Initialize counter*/
      int count=0;
 
      /*Create spcs from x,y,z, as well as the spc values on those spcs: */
      for(int i=0;i<iomodel->numberofvertices;i++){
 
         /*keep only this partition's nodes:*/
         if(iomodel->my_vertices[i]){
 
            if (reCast<int,IssmDouble>(iomodel->Data("md.flowequation.borderFS")[i])){
               constraints->AddObject(new SpcStatic(count+1,i+1,0,0,StressbalanceVerticalAnalysisEnum)); //spc to zero as vertical velocity is done in Horiz for FS
               count++;
            }
            else if (!xIsNan<IssmDouble>(spcvz[i])){
               constraints->AddObject(new SpcStatic(count+1,i+1,0,
                           spcvz[i],StressbalanceVerticalAnalysisEnum)); //add count'th spc, on node i+1, setting dof 1 to vx.
               count++;
 
            }
         } 
      }
 
      /*Free data: */
      iomodel->DeleteData(1,"md.flowequation.borderFS");
      iomodel->DeleteData(spcvz,"md.stressbalance.spcvz");
   }
 
}/*}}}*/

CreateLoads()

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

virtual

Implements Analysis.

Definition at line 71 of file StressbalanceVerticalAnalysis.cpp.

                                                                             {/*{{{*/
 
   /*No loads*/
 
}/*}}}*/

CreateNodes()

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

virtual

Implements Analysis.

Definition at line 76 of file StressbalanceVerticalAnalysis.cpp.

                                                                                       {/*{{{*/
 
   /*return if not 3d mesh*/
   if(iomodel->domaintype!=Domain3DEnum) return;
 
   iomodel->FetchData(3,"md.mesh.vertexonbase","md.mesh.vertexonsurface","md.flowequation.vertex_equation");
   ::CreateNodes(nodes,iomodel,StressbalanceVerticalAnalysisEnum,P1Enum);
   iomodel->DeleteData(3,"md.mesh.vertexonbase","md.mesh.vertexonsurface","md.flowequation.vertex_equation");
}/*}}}*/

DofsPerNode()

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

virtual

Implements Analysis.

Definition at line 85 of file StressbalanceVerticalAnalysis.cpp.

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

UpdateElements()

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

virtual

Implements Analysis.

Definition at line 88 of file StressbalanceVerticalAnalysis.cpp.

                                                                                                                                             {/*{{{*/
 
   /*return if not 3d mesh*/
   if(iomodel->domaintype!=Domain3DEnum) return;
 
   /*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,P1Enum);
         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);
   if(iomodel->domaintype!=Domain2DhorizontalEnum){
      iomodel->FetchDataToInput(inputs2,elements,"md.mesh.vertexonbase",MeshVertexonbaseEnum);
      iomodel->FetchDataToInput(inputs2,elements,"md.mesh.vertexonsurface",MeshVertexonsurfaceEnum);
   }
   iomodel->FetchDataToInput(inputs2,elements,"md.basalforcings.groundedice_melting_rate",BasalforcingsGroundediceMeltingRateEnum);
   //iomodel->FetchDataToInput(inputs2,elements,"md.smb.mass_balance",SmbMassBalanceEnum);
 
 
   /*Add basal forcings to compute melt rate*/
   int basalforcing_model;
   iomodel->FindConstant(&basalforcing_model,"md.basalforcings.model");
   switch(basalforcing_model){
      case FloatingMeltRateEnum:
         iomodel->FetchDataToInput(inputs2,elements,"md.basalforcings.floatingice_melting_rate",BasalforcingsFloatingiceMeltingRateEnum);
         break;
      case LinearFloatingMeltRateEnum:
         break;
      case MismipFloatingMeltRateEnum:
         break;
      case MantlePlumeGeothermalFluxEnum:
         break;
      case SpatialLinearFloatingMeltRateEnum:
         iomodel->FetchDataToInput(inputs2,elements,"md.basalforcings.deepwater_melting_rate",BasalforcingsDeepwaterMeltingRateEnum);
         iomodel->FetchDataToInput(inputs2,elements,"md.basalforcings.deepwater_elevation",BasalforcingsDeepwaterElevationEnum);
         iomodel->FetchDataToInput(inputs2,elements,"md.basalforcings.upperwater_elevation",BasalforcingsUpperwaterElevationEnum);
         break;
      case BasalforcingsPicoEnum:
         iomodel->FetchDataToInput(inputs2,elements,"md.basalforcings.basin_id",BasalforcingsPicoBasinIdEnum);
         break;
      case BasalforcingsIsmip6Enum:
         iomodel->FetchDataToInput(inputs2,elements,"md.basalforcings.basin_id",BasalforcingsIsmip6BasinIdEnum);
         break;
      case BeckmannGoosseFloatingMeltRateEnum:
         iomodel->FetchDataToInput(inputs2,elements,"md.basalforcings.ocean_salinity",BasalforcingsOceanSalinityEnum);
         iomodel->FetchDataToInput(inputs2,elements,"md.basalforcings.ocean_temp",BasalforcingsOceanTempEnum);
         break;
      default:
         _error_("Basal forcing model "<<EnumToStringx(basalforcing_model)<<" not supported yet");
   }
   iomodel->FetchDataToInput(inputs2,elements,"md.initialization.vx",VxEnum,0.);
   iomodel->FetchDataToInput(inputs2,elements,"md.initialization.vy",VyEnum,0.);
}/*}}}*/

UpdateParameters()

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

virtual

Implements Analysis.

Definition at line 150 of file StressbalanceVerticalAnalysis.cpp.

                                                                                                                               {/*{{{*/
 
   /*No specific parameters*/
 
}/*}}}*/

Core()

void StressbalanceVerticalAnalysis::Core ( FemModel *  femmodel )

virtual

Implements Analysis.

Definition at line 157 of file StressbalanceVerticalAnalysis.cpp.

                                                                    {/*{{{*/
 
      if(VerboseSolution()) _printf0_("   computing vertical velocities\n");
      femmodel->SetCurrentConfiguration(StressbalanceVerticalAnalysisEnum);
      solutionsequence_linear(femmodel);
}/*}}}*/

CreateDVector()

ElementVector * StressbalanceVerticalAnalysis::CreateDVector ( Element *  element )

virtual

Implements Analysis.

Definition at line 163 of file StressbalanceVerticalAnalysis.cpp.

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

CreateJacobianMatrix()

ElementMatrix * StressbalanceVerticalAnalysis::CreateJacobianMatrix ( Element *  element )

virtual

Implements Analysis.

Definition at line 167 of file StressbalanceVerticalAnalysis.cpp.

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

CreateKMatrix()

ElementMatrix * StressbalanceVerticalAnalysis::CreateKMatrix ( Element *  element )

virtual

Implements Analysis.

Definition at line 170 of file StressbalanceVerticalAnalysis.cpp.

                                                                           {/*{{{*/
 
   /* Check if ice in element */
   if(!element->IsIceInElement()) return NULL;
 
   bool hack = false;
 
   /*compute all stiffness matrices for this element*/
   ElementMatrix* Ke1=CreateKMatrixVolume(element);
   ElementMatrix* Ke2=NULL;
   if(hack) Ke2=CreateKMatrixBase(element);
   else Ke2=CreateKMatrixSurface(element);
   ElementMatrix* Ke =new ElementMatrix(Ke1,Ke2);
 
   /*clean-up and return*/
   delete Ke1;
   delete Ke2;
   return Ke;
 
}/*}}}*/

CreateKMatrixBase()

ElementMatrix * StressbalanceVerticalAnalysis::CreateKMatrixBase

(

Element *

element

)

Definition at line 190 of file StressbalanceVerticalAnalysis.cpp.

                                                                               {/*{{{*/
 
   if(!element->IsOnBase()) return NULL;
 
   /*Intermediaries*/
   IssmDouble  D,Jdet,normal[3];
   IssmDouble *xyz_list = NULL;
 
   /*Fetch number of nodes and dof for this finite element*/
   int numnodes = element->GetNumberOfNodes();
 
   /*Initialize Element matrix and vectors*/
   ElementMatrix* Ke    = element->NewElementMatrix(NoneApproximationEnum);
   IssmDouble*    basis = xNew<IssmDouble>(numnodes);
 
   /*Retrieve all inputs and parameters*/
   element->GetVerticesCoordinatesBase(&xyz_list);
 
   /* Start  looping on the number of gaussian points: */
   Gauss* gauss = element->NewGaussBase(2);
   element->NormalBase(&normal[0],xyz_list);
   for(int ig=gauss->begin();ig<gauss->end();ig++){
      gauss->GaussPoint(ig);
 
      element->JacobianDeterminantBase(&Jdet,xyz_list,gauss);
      element->NodalFunctions(basis,gauss);
      D = -gauss->weight*Jdet*normal[2];
 
      for(int i=0;i<numnodes;i++) for(int j=0;j<numnodes;j++) Ke->values[i*numnodes+j] += D*basis[i]*basis[j];
   }
 
   /*Clean up and return*/
   delete gauss;
   xDelete<IssmDouble>(xyz_list);
   xDelete<IssmDouble>(basis);
   return Ke;
}/*}}}*/

CreateKMatrixSurface()

ElementMatrix * StressbalanceVerticalAnalysis::CreateKMatrixSurface

(

Element *

element

)

Definition at line 227 of file StressbalanceVerticalAnalysis.cpp.

                                                                                  {/*{{{*/
 
   if(!element->IsOnSurface()) return NULL;
 
   /*Intermediaries*/
   IssmDouble  D,Jdet,normal[3];
   IssmDouble *xyz_list = NULL;
 
   /*Fetch number of nodes and dof for this finite element*/
   int numnodes = element->GetNumberOfNodes();
 
   /*Initialize Element matrix and vectors*/
   ElementMatrix* Ke    = element->NewElementMatrix(NoneApproximationEnum);
   IssmDouble*    basis = xNew<IssmDouble>(numnodes);
 
   /*Retrieve all inputs and parameters*/
   element->GetVerticesCoordinatesTop(&xyz_list);
 
   /* Start  looping on the number of gaussian points: */
   Gauss* gauss = element->NewGaussTop(2);
   element->NormalTop(&normal[0],xyz_list);
   for(int ig=gauss->begin();ig<gauss->end();ig++){
      gauss->GaussPoint(ig);
 
      element->JacobianDeterminantTop(&Jdet,xyz_list,gauss);
      element->NodalFunctions(basis,gauss);
      D = -gauss->weight*Jdet*normal[2];
 
      for(int i=0;i<numnodes;i++) for(int j=0;j<numnodes;j++) Ke->values[i*numnodes+j] += D*basis[i]*basis[j];
   }
 
   /*Clean up and return*/
   delete gauss;
   xDelete<IssmDouble>(xyz_list);
   xDelete<IssmDouble>(basis);
   return Ke;
}/*}}}*/

CreateKMatrixVolume()

ElementMatrix * StressbalanceVerticalAnalysis::CreateKMatrixVolume

(

Element *

element

)

Definition at line 264 of file StressbalanceVerticalAnalysis.cpp.

                                                                                 {/*{{{*/
 
   /*Intermediaries*/
   IssmDouble  Jdet;
   IssmDouble *xyz_list = NULL;
 
   /*Fetch number of nodes and dof for this finite element*/
   int numnodes = element->GetNumberOfNodes();
 
   /*Initialize Element matrix and vectors*/
   ElementMatrix* Ke     = element->NewElementMatrix(NoneApproximationEnum);
   IssmDouble*    basis  = xNew<IssmDouble>(numnodes);
   IssmDouble*    dbasis = xNew<IssmDouble>(3*numnodes);
 
   /*Retrieve all inputs and parameters*/
   element->GetVerticesCoordinates(&xyz_list);
 
   /* 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);
 
      for(int i=0;i<numnodes;i++){
         for(int j=0;j<numnodes;j++){
            Ke->values[i*numnodes+j] += gauss->weight*Jdet*(
                     basis[j]*dbasis[2*numnodes+i]
                     );
         }
      }
   }
 
   /*Clean up and return*/
   delete gauss;
   xDelete<IssmDouble>(xyz_list);
   xDelete<IssmDouble>(dbasis);
   xDelete<IssmDouble>(basis);
   return Ke;
 
}/*}}}*/

CreatePVector()

ElementVector * StressbalanceVerticalAnalysis::CreatePVector ( Element *  element )

virtual

Implements Analysis.

Definition at line 307 of file StressbalanceVerticalAnalysis.cpp.

                                                                           {/*{{{*/
 
   /* Check if ice in element */
   if(!element->IsIceInElement()) return NULL;
 
   bool hack = false;
 
   /*compute all load vectors for this element*/
   ElementVector* pe1=CreatePVectorVolume(element);
   ElementVector* pe2=NULL;
   if(hack) pe2=CreatePVectorSurface(element);
   else     pe2=CreatePVectorBase(element);
   ElementVector* pe =new ElementVector(pe1,pe2);
 
   /*clean-up and return*/
   delete pe1;
   delete pe2;
   return pe;
}/*}}}*/

CreatePVectorBase()

ElementVector * StressbalanceVerticalAnalysis::CreatePVectorBase

(

Element *

element

)

Definition at line 326 of file StressbalanceVerticalAnalysis.cpp.

                                                                               {/*{{{*/
 
   /*Intermediaries */
   int         approximation;
   IssmDouble *xyz_list      = NULL;
   IssmDouble *xyz_list_base = NULL;
   IssmDouble  Jdet,slope[3];
   IssmDouble  vx,vy,vz=0.,dbdx,dbdy;
   IssmDouble  gmb,fmb,phi,basalmeltingvalue;
 
   if(!element->IsOnBase()) return NULL;
 
   /*Fetch number of nodes for this finite element*/
   int numnodes = element->GetNumberOfNodes();
 
   /*Initialize Element vector*/
   ElementVector* pe    = element->NewElementVector();
   IssmDouble*    basis = xNew<IssmDouble>(numnodes);
 
   /*Retrieve all inputs and parameters*/
   element->GetVerticesCoordinates(&xyz_list);
   element->GetVerticesCoordinatesBase(&xyz_list_base);
   element->GetInputValue(&approximation,ApproximationEnum);
   Input2* base_input=element->GetInput2(BaseEnum);                                               _assert_(base_input);
   Input2* groundedice_input=element->GetInput2(MaskOceanLevelsetEnum);                     _assert_(groundedice_input);
   Input2* groundedice_melting_input=element->GetInput2(BasalforcingsGroundediceMeltingRateEnum); _assert_(groundedice_melting_input);
   Input2* floatingice_melting_input=element->GetInput2(BasalforcingsFloatingiceMeltingRateEnum); _assert_(floatingice_melting_input);
   Input2* vx_input=element->GetInput2(VxEnum);                                                   _assert_(vx_input);
   Input2* vy_input=element->GetInput2(VyEnum);                                                   _assert_(vy_input);
   Input2* vzFS_input=NULL;
   if(approximation==HOFSApproximationEnum || approximation==SSAFSApproximationEnum){
      vzFS_input=element->GetInput2(VzFSEnum);       _assert_(vzFS_input);
   }
 
   /* Start  looping on the number of gaussian points: */
   Gauss* gauss=element->NewGaussBase(2);
   for(int ig=gauss->begin();ig<gauss->end();ig++){
      gauss->GaussPoint(ig);
 
      groundedice_melting_input->GetInputValue(&gmb,gauss);
      floatingice_melting_input->GetInputValue(&fmb,gauss);
      groundedice_input->GetInputValue(&phi,gauss);
      base_input->GetInputDerivativeValue(&slope[0],xyz_list,gauss);
      vx_input->GetInputValue(&vx,gauss);
      vy_input->GetInputValue(&vy,gauss);
      if(approximation==HOFSApproximationEnum || approximation==SSAFSApproximationEnum){
         vzFS_input->GetInputValue(&vz,gauss);
      }
      dbdx=slope[0];
      dbdy=slope[1];
      if(phi>0.) basalmeltingvalue=gmb;
      else basalmeltingvalue=fmb;
 
      element->JacobianDeterminantBase(&Jdet,xyz_list_base,gauss);
      element->NodalFunctions(basis,gauss);
 
      for(int i=0;i<numnodes;i++) pe->values[i]+=-Jdet*gauss->weight*(vx*dbdx+vy*dbdy-vz-basalmeltingvalue)*basis[i];
   }
 
   /*Clean up and return*/
   delete gauss;
   xDelete<IssmDouble>(basis);
   xDelete<IssmDouble>(xyz_list);
   xDelete<IssmDouble>(xyz_list_base);
   return pe;
}/*}}}*/

CreatePVectorSurface()

ElementVector * StressbalanceVerticalAnalysis::CreatePVectorSurface

(

Element *

element

)

Definition at line 392 of file StressbalanceVerticalAnalysis.cpp.

                                                                                  {/*{{{*/
 
   /*Intermediaries */
   int         approximation;
   IssmDouble *xyz_list      = NULL;
   IssmDouble *xyz_list_surface= NULL;
   IssmDouble  Jdet,slope[3];
   IssmDouble  vx,vy,vz=0.,dsdx,dsdy;
   IssmDouble  smb,smbvalue;
 
   if(!element->IsOnSurface()) return NULL;
 
   /*Fetch number of nodes for this finite element*/
   int numnodes = element->GetNumberOfNodes();
 
   /*Initialize Element vector*/
   ElementVector* pe    = element->NewElementVector();
   IssmDouble*    basis = xNew<IssmDouble>(numnodes);
 
   /*Retrieve all inputs and parameters*/
   element->GetVerticesCoordinates(&xyz_list);
   element->GetVerticesCoordinatesTop(&xyz_list_surface);
   element->GetInputValue(&approximation,ApproximationEnum);
   Input2* surface_input    =element->GetInput2(SurfaceEnum);   _assert_(surface_input);
   Input2* smb_input=element->GetInput2(SmbMassBalanceEnum);    _assert_(smb_input);
   Input2* vx_input=element->GetInput2(VxEnum);                 _assert_(vx_input);
   Input2* vy_input=element->GetInput2(VyEnum);                 _assert_(vy_input);
   Input2* vzFS_input=NULL;
   if(approximation==HOFSApproximationEnum || approximation==SSAFSApproximationEnum){
      vzFS_input=element->GetInput2(VzFSEnum); _assert_(vzFS_input);
   }
 
   /* Start  looping on the number of gaussian points: */
   Gauss* gauss=element->NewGaussTop(2);
   for(int ig=gauss->begin();ig<gauss->end();ig++){
      gauss->GaussPoint(ig);
 
      smb_input->GetInputValue(&smb,gauss);
      surface_input->GetInputDerivativeValue(&slope[0],xyz_list,gauss);
      vx_input->GetInputValue(&vx,gauss);
      vy_input->GetInputValue(&vy,gauss);
      if(approximation==HOFSApproximationEnum || approximation==SSAFSApproximationEnum){
         vzFS_input->GetInputValue(&vz,gauss);
      }
      dsdx=slope[0];
      dsdy=slope[1];
 
      element->JacobianDeterminantTop(&Jdet,xyz_list_surface,gauss);
      element->NodalFunctions(basis,gauss);
 
      for(int i=0;i<numnodes;i++) pe->values[i]+=-Jdet*gauss->weight*(vx*dsdx+vy*dsdy-vz+smb)*basis[i];
   }
 
   /*Clean up and return*/
   delete gauss;
   xDelete<IssmDouble>(basis);
   xDelete<IssmDouble>(xyz_list);
   xDelete<IssmDouble>(xyz_list_surface);
   return pe;
}/*}}}*/

CreatePVectorVolume()

ElementVector * StressbalanceVerticalAnalysis::CreatePVectorVolume

(

Element *

element

)

Definition at line 452 of file StressbalanceVerticalAnalysis.cpp.

                                                                                 {/*{{{*/
 
   /*Intermediaries*/
   int         approximation;
   IssmDouble  Jdet,dudx,dvdy,dwdz;
   IssmDouble  du[3],dv[3],dw[3];
   IssmDouble* xyz_list = NULL;
 
   /*Fetch number of nodes for this finite element*/
   int numnodes = element->GetNumberOfNodes();
 
   /*Initialize Element vector and basis functions*/
   ElementVector* pe    = element->NewElementVector();
   IssmDouble*    basis = xNew<IssmDouble>(numnodes);
 
   /*Retrieve all inputs and parameters*/
   element->GetVerticesCoordinates(&xyz_list);
   element->GetInputValue(&approximation,ApproximationEnum);
   Input2* vx_input=element->GetInput2(VxEnum); _assert_(vx_input);
   Input2* vy_input=element->GetInput2(VyEnum); _assert_(vy_input);
   Input2* vzFS_input=NULL;
   if(approximation==HOFSApproximationEnum || approximation==SSAFSApproximationEnum){
      vzFS_input=element->GetInput2(VzFSEnum); _assert_(vzFS_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);
      element->NodalFunctions(basis,gauss);
 
      vx_input->GetInputDerivativeValue(&du[0],xyz_list,gauss);
      vy_input->GetInputDerivativeValue(&dv[0],xyz_list,gauss);
      if(approximation==HOFSApproximationEnum || approximation==SSAFSApproximationEnum){
         vzFS_input->GetInputDerivativeValue(&dw[0],xyz_list,gauss);
         dwdz=dw[2];
      }
      else dwdz=0;
      dudx=du[0];
      dvdy=dv[1];
 
      for(int i=0;i<numnodes;i++) pe->values[i] += (dudx+dvdy+dwdz)*Jdet*gauss->weight*basis[i];
   }
 
   /*Clean up and return*/
   delete gauss;
   xDelete<IssmDouble>(basis);
   xDelete<IssmDouble>(xyz_list);
   return pe;
}/*}}}*/

GetSolutionFromInputs()

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

virtual

Implements Analysis.

Definition at line 504 of file StressbalanceVerticalAnalysis.cpp.

                                                                                                                {/*{{{*/
   element->GetSolutionFromInputsOneDof(solution,VzEnum);
}/*}}}*/

GradientJ()

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

virtual

Implements Analysis.

Definition at line 507 of file StressbalanceVerticalAnalysis.cpp.

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

InputUpdateFromSolution()

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

virtual

Implements Analysis.

Definition at line 510 of file StressbalanceVerticalAnalysis.cpp.

                                                                                                          {/*{{{*/
 
   int          numnodes = element->GetNumberOfNodes();
   int          numdof=numnodes*1;
 
   int          i;
   int          approximation;
   int*         doflist  = NULL;
   IssmDouble*  xyz_list = NULL;
   IssmDouble   rho_ice,g;
 
   /*Get the approximation and do nothing if the element in FS or None*/
   element->GetInputValue(&approximation,ApproximationEnum);
   if(approximation==FSApproximationEnum || approximation==NoneApproximationEnum){
      return;
   }
 
   /*Get dof list and vertices coordinates: */
   element->GetVerticesCoordinates(&xyz_list);
   element->GetDofListLocal(&doflist,NoneApproximationEnum,GsetEnum);
   IssmDouble*  values    = xNew<IssmDouble>(numdof);
   IssmDouble*  vx        = xNew<IssmDouble>(numnodes);
   IssmDouble*  vy        = xNew<IssmDouble>(numnodes);
   IssmDouble*  vz        = xNew<IssmDouble>(numnodes);
   IssmDouble*  vzSSA     = xNew<IssmDouble>(numnodes);
   IssmDouble*  vzHO      = xNew<IssmDouble>(numnodes);
   IssmDouble*  vzFS      = xNew<IssmDouble>(numnodes);
   IssmDouble*  vel       = xNew<IssmDouble>(numnodes);
   IssmDouble*  pressure  = xNew<IssmDouble>(numnodes);
   IssmDouble*  surface   = xNew<IssmDouble>(numnodes);
 
   /*Use the dof list to index into the solution vector vz: */
   for(i=0;i<numdof;i++) values[i]=solution[doflist[i]];
   for(i=0;i<numdof;i++){
      vz[i]=values[i*1+0];
 
      /*Check solution*/
      if(xIsNan<IssmDouble>(vz[i])) _error_("NaN found in solution vector");
      if(xIsInf<IssmDouble>(vz[i])) _error_("Inf found in solution vector");
   }
 
   /*Get Vx and Vy*/
   element->GetInputListOnNodes(&vx[0],VxEnum,0.0); //default is 0
   element->GetInputListOnNodes(&vy[0],VyEnum,0.0); //default is 0
 
   /*Do some modifications if we actually have a HOFS or SSAFS element*/
   if(approximation==HOFSApproximationEnum){
      Input2* vzFS_input=element->GetInput2(VzFSEnum);
      if (vzFS_input){
         if (vzFS_input->ObjectEnum()!=PentaInput2Enum) _error_("Cannot compute Vel as VzFS is of type " << EnumToStringx(vzFS_input->ObjectEnum()));
         element->GetInputListOnNodes(&vzFS[0],VzFSEnum,0.);
      }
      else _error_("Cannot compute Vz as VzFS in not present in HOFS element");
      for(i=0;i<numnodes;i++){
         vzHO[i]=vz[i];
         vz[i]=vzHO[i]+vzFS[i];
      }
   }
   else if(approximation==SSAFSApproximationEnum){
      Input2* vzFS_input=element->GetInput2(VzFSEnum);
      if (vzFS_input){
         if (vzFS_input->ObjectEnum()!=PentaInput2Enum) _error_("Cannot compute Vel as VzFS is of type " << EnumToStringx(vzFS_input->ObjectEnum()));
         element->GetInputListOnNodes(&vzFS[0],VzFSEnum,0.);
      }
      else _error_("Cannot compute Vz as VzFS in not present in SSAFS element");
      for(i=0;i<numnodes;i++){
         vzSSA[i]=vz[i];
         vz[i]=vzSSA[i]+vzFS[i];
      }
   }
 
   /*Now Compute vel*/
   for(i=0;i<numnodes;i++) vel[i]=pow( pow(vx[i],2.0) + pow(vy[i],2.0) + pow(vz[i],2.0) , 0.5);
 
   /*For pressure: we have not computed pressure in this analysis, for this element. We are in 3D, 
    *so the pressure is just the pressure at the z elevation: except it this is a HOFS element */
   if(approximation!=HOFSApproximationEnum &&  approximation!=SSAFSApproximationEnum){
      rho_ice = element->FindParam(MaterialsRhoIceEnum);
      g       = element->FindParam(ConstantsGEnum);
      element->GetInputListOnNodes(&surface[0],SurfaceEnum,0.);
      for(i=0;i<numnodes;i++) pressure[i]=rho_ice*g*(surface[i]-xyz_list[i*3+2]);
   }
   if(approximation!=HOFSApproximationEnum && approximation!=SSAFSApproximationEnum){
      element->AddInput2(PressureEnum,pressure,element->GetElementType());
   }
   else if(approximation==HOFSApproximationEnum){
      element->AddInput2(VzHOEnum,vzHO,P1Enum);
   }
   else if(approximation==SSAFSApproximationEnum){
      element->AddInput2(VzSSAEnum,vzSSA,P1Enum);
   }
   element->AddInput2(VzEnum,vz,P1Enum);
   element->AddInput2(VelEnum,vel,P1Enum);
 
   /*Free ressources:*/
   xDelete<IssmDouble>(surface);
   xDelete<IssmDouble>(pressure);
   xDelete<IssmDouble>(vel);
   xDelete<IssmDouble>(vz);
   xDelete<IssmDouble>(vzSSA);
   xDelete<IssmDouble>(vzHO);
   xDelete<IssmDouble>(vzFS);
   xDelete<IssmDouble>(vy);
   xDelete<IssmDouble>(vx);
   xDelete<IssmDouble>(values);
   xDelete<IssmDouble>(xyz_list);
   xDelete<int>(doflist);
}/*}}}*/

UpdateConstraints()

void StressbalanceVerticalAnalysis::UpdateConstraints ( FemModel *  femmodel )

virtual

Implements Analysis.

Definition at line 618 of file StressbalanceVerticalAnalysis.cpp.

                                                                                 {/*{{{*/
   /*Default, do nothing*/
   SetActiveNodesLSMx(femmodel);
   return;
}/*}}}*/

---

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

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