Cfsurfacesquare.cpp

Go to the documentation of this file.

 
/*Headers:*/
/*{{{*/
#ifdef HAVE_CONFIG_H
   #include <config.h>
#else
#error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!"
#endif
 
#include "./classes.h"
#include "./ExternalResults/ExternalResult.h"
#include "./ExternalResults/Results.h"
#include "../datastructures/datastructures.h"
#include "./Elements/Element.h"
#include "./Elements/Elements.h"
#include "./FemModel.h"
#include "../modules/SurfaceAreax/SurfaceAreax.h"
#include "../classes/Params/Parameters.h"
#include "../classes/gauss/Gauss.h"
#include "./Inputs2/DatasetInput2.h"
/*}}}*/
 
/*Cfsurfacesquare constructors, destructors :*/
Cfsurfacesquare::Cfsurfacesquare(){/*{{{*/
 
   this->definitionenum = -1;
   this->name = NULL;
   this->model_enum = UNDEF;
   this->observation_enum = UNDEF;
   this->weights_enum = UNDEF;
   this->misfit=0;
   this->lock=0;
   this->datatime=0.;
   this->timepassedflag = false;
   this->last_time = 0.;
 
}
/*}}}*/
Cfsurfacesquare::Cfsurfacesquare(char* in_name, int in_definitionenum, int in_model_enum, int in_observation_enum, int in_weights_enum, IssmDouble in_datatime, bool in_timepassedflag){/*{{{*/
 
   this->definitionenum=in_definitionenum;
 
   this->name     = xNew<char>(strlen(in_name)+1);
   xMemCpy<char>(this->name,in_name,strlen(in_name)+1);
 
   this->model_enum=in_model_enum;
   this->observation_enum=in_observation_enum;
   this->weights_enum=in_weights_enum;
   this->datatime=in_datatime;
   this->timepassedflag=in_timepassedflag;
 
   this->misfit=0;
   this->lock=0;
}
/*}}}*/
Cfsurfacesquare::~Cfsurfacesquare(){/*{{{*/
   if(this->name)xDelete(this->name);
   this->misfit=0;
   this->lock=0;
}
/*}}}*/
/*Object virtual function resolutoin: */
Object* Cfsurfacesquare::copy() {/*{{{*/
   Cfsurfacesquare* mf = new Cfsurfacesquare(this->name,this->definitionenum, this->model_enum,this->observation_enum,this->weights_enum,this->datatime,this->timepassedflag);
   mf->misfit=this->misfit;
   mf->lock=this->lock;
   return (Object*) mf;
}
/*}}}*/
void Cfsurfacesquare::DeepEcho(void){/*{{{*/
   this->Echo();
}
/*}}}*/
void Cfsurfacesquare::Echo(void){/*{{{*/
   _printf_(" Cfsurfacesquare: " << name << " " << this->definitionenum << "\n");
   _printf_("    model_enum: " << model_enum << " " << EnumToStringx(model_enum) << "\n");
   _printf_("    observation_enum: " << observation_enum << " " << EnumToStringx(observation_enum) << "\n");
   _printf_("    weights_enum: " << weights_enum << " " << EnumToStringx(weights_enum) << "\n");
   _printf_("    datatime: " << datatime << "\n");
   _printf_("    timepassedflag: "<<timepassedflag<<"\n");
}
/*}}}*/
int Cfsurfacesquare::Id(void){/*{{{*/
   return -1;
}
/*}}}*/
void Cfsurfacesquare::Marshall(char** pmarshalled_data,int* pmarshalled_data_size, int marshall_direction){/*{{{*/
   _error_("not implemented yet!"); 
} 
/*}}}*/
int Cfsurfacesquare::ObjectEnum(void){/*{{{*/
   return CfsurfacesquareEnum;
}
/*}}}*/
/*Definition virtual function resolutoin: */
int Cfsurfacesquare::DefinitionEnum(){/*{{{*/
   return this->definitionenum;
}
/*}}}*/
char* Cfsurfacesquare::Name(){/*{{{*/
   char* name2=xNew<char>(strlen(this->name)+1);
   xMemCpy(name2,this->name,strlen(this->name)+1);
 
   return name2;
}
/*}}}*/
IssmDouble Cfsurfacesquare::Response(FemModel* femmodel){/*{{{*/
    /*diverse: */
    IssmDouble time;
 
    /*recover time parameters: */
    femmodel->parameters->FindParam(&time,TimeEnum);
    if(time < last_time) timepassedflag = false;
    last_time = time;
 
       int i;
       IssmDouble J=0.;
       IssmDouble J_sum=0.;
 
    if(datatime<=time && !timepassedflag){
       for(i=0;i<femmodel->elements->Size();i++){
          Element* element=(Element*)femmodel->elements->GetObjectByOffset(i);
          J+=this->Cfsurfacesquare_Calculation(element,model_enum,observation_enum,weights_enum);
       }
 
       ISSM_MPI_Allreduce ( (void*)&J,(void*)&J_sum,1,ISSM_MPI_DOUBLE,ISSM_MPI_SUM,IssmComm::GetComm());
       ISSM_MPI_Bcast(&J_sum,1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm());
       J=J_sum;
 
       timepassedflag = true;
       return J;
   }
    else return J;
 }
   /*}}}*/
IssmDouble Cfsurfacesquare::Cfsurfacesquare_Calculation(Element* element, int model_enum, int observation_enum, int weights_enum){/*{{{*/
 
   int        domaintype,numcomponents;
   IssmDouble Jelem=0.;
   IssmDouble misfit,Jdet;
   IssmDouble model,obs,weight;
   IssmDouble* xyz_list = NULL;
 
   /*Get basal element*/
   if(!element->IsOnSurface()) return 0.;
 
   /*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");
   }
 
   /*Spawn surface element*/
   Element* topelement = element->SpawnTopElement();
 
   /* Get node coordinates*/
   topelement->GetVerticesCoordinates(&xyz_list);
 
   /*Retrieve all inputs we will be needing: */
   DatasetInput2 *datasetinput = topelement->GetDatasetInput2(definitionenum); _assert_(datasetinput);
   Input2        *model_input  = topelement->GetInput2(model_enum);            _assert_(model_input);
 
   /* Start  looping on the number of gaussian points: */
   Gauss* gauss=topelement->NewGauss(2);
   for(int ig=gauss->begin();ig<gauss->end();ig++){
 
      gauss->GaussPoint(ig);
 
      /* Get Jacobian determinant: */
      topelement->JacobianDeterminant(&Jdet,xyz_list,gauss);
 
      /*Get all parameters at gaussian point*/
      datasetinput->GetInputValue(&weight,gauss,WeightsSurfaceObservationEnum);
      model_input->GetInputValue(&model,gauss);
      datasetinput->GetInputValue(&obs,gauss,SurfaceObservationEnum);
 
      /*Compute SurfaceAbsVelMisfitEnum:
       *        *
       *               *      1  [           2              2 ]
       *                      * J = --- | (u - u   )  +  (v - v   )  |
       *                             *      2  [       obs            obs   ]
       *                                    *
       *                                           */
      misfit=0.5*(model-obs)*(model-obs);
 
      /*Add to cost function*/
      Jelem+=misfit*weight*Jdet*gauss->weight;
   }
 
   /*clean up and Return: */
   if(domaintype!=Domain2DhorizontalEnum){topelement->DeleteMaterials(); delete topelement;};
   xDelete<IssmDouble>(xyz_list);
   delete gauss;
   return Jelem;
}/*}}}*/