gia_core.cpp

Go to the documentation of this file.

 
#include "./cores.h"
#include "../toolkits/toolkits.h"
#include "../classes/classes.h"
#include "../classes/Inputs2/TriaInput2.h"
#include "../classes/Inputs2/DatasetInput2.h"
#include "../shared/shared.h"
#include "../modules/modules.h"
#include "../solutionsequences/solutionsequences.h"
void gia_core(FemModel* femmodel){
 
   Vector<IssmDouble> *wg    = NULL;
   Vector<IssmDouble> *dwdtg = NULL;
   IssmDouble          *x    = NULL;
   IssmDouble          *y    = NULL;
   int nummodels,giamodel;
   IssmDouble modelid; 
 
   /*parameters: */
   bool save_results;
   int  gsize;
   
   /*Start profiler*/
   femmodel->profiler->Start(GIACORE);
 
   /*Recover some parameters: */
   femmodel->parameters->FindParam(&save_results,SaveResultsEnum);
 
   if(VerboseSolution()) _printf0_("   computing GIA\n");
 
   /*Call on core computations: */
   femmodel->SetCurrentConfiguration(GiaAnalysisEnum);
 
   /*Figure out type of gia model:*/
   femmodel->parameters->FindParam(&giamodel,GiaModelEnum);
 
   if(giamodel==1){ //GiaIvins
      /*Figure out size of g-set deflection vector and allocate solution vector: */
      gsize      = femmodel->nodes->NumberOfDofs(GsetEnum);
      wg = new Vector<IssmDouble>(gsize);
      dwdtg = new Vector<IssmDouble>(gsize);
 
      /*first, recover x and y vectors from vertices: */
      VertexCoordinatesx(&x,&y,NULL,femmodel->vertices); //no need for z coordinate
 
      /*call the main module: */
      femmodel->Deflection(wg,dwdtg,x,y);
 
      /*assemble vector: */
      wg->Assemble();
      dwdtg->Assemble();
 
      InputUpdateFromVectorx(femmodel,wg,UGiaEnum,VertexSIdEnum);
      InputUpdateFromVectorx(femmodel,dwdtg,UGiaRateEnum,VertexSIdEnum);
 
      if(save_results){
         if(VerboseSolution()) _printf0_("   saving results\n");
         int outputs[2] = {UGiaEnum,UGiaRateEnum};
         femmodel->RequestedOutputsx(&femmodel->results,&outputs[0],2);
      }
 
      xDelete<IssmDouble>(x);
      xDelete<IssmDouble>(y);
   }
   else if(giamodel==2){ //GiaCaron
      /*not implemneted yet*/
   }
   else if(giamodel==3){ //GiaMme
      
      /*Go grab the values for Ngia and Ugia in the offline GiaMme class:*/ 
      femmodel->parameters->FindParam(&modelid,GiaModelidEnum); 
      modelid--; //from matlab
 
      /*find the Ngia and Ugia dataset:*/
      DatasetInput2* dataset_input_ngia=femmodel->inputs2->GetDatasetInput2(GiaMmeNgiaEnum);
      DatasetInput2* dataset_input_ugia=femmodel->inputs2->GetDatasetInput2(GiaMmeUgiaEnum);
      /*Go find the modelid'th input:*/
      TriaInput2* tria_input_ngia=dataset_input_ngia->GetTriaInputByOffset(modelid);
      TriaInput2* tria_input_ugia=dataset_input_ugia->GetTriaInputByOffset(modelid);
 
      /*Plug into SealevelUGiaRate and SealevelNGiaRate inputs:*/
      Input2* tria_input_copy_ngia=tria_input_ngia->copy();
      Input2* tria_input_copy_ugia=tria_input_ugia->copy();
      tria_input_copy_ngia->ChangeEnum(NGiaRateEnum);
      tria_input_copy_ugia->ChangeEnum(UGiaRateEnum);
      femmodel->inputs2->AddInput(tria_input_copy_ngia);
      femmodel->inputs2->AddInput(tria_input_copy_ugia);
   }
   
   /*End profiler*/
   femmodel->profiler->Stop(GIACORE);
}