Index: /issm/trunk/src/c/Container/Inputs.cpp =================================================================== --- /issm/trunk/src/c/Container/Inputs.cpp (revision 8362) +++ /issm/trunk/src/c/Container/Inputs.cpp (revision 8363) @@ -429,2 +429,17 @@ } /*}}}*/ +/*FUNCTION Inputs::Configure{{{1*/ +void Inputs::Configure(Parameters* parameters){ + + vector::iterator object; + Input* input=NULL; + + for ( object=objects.begin() ; object < objects.end(); object++ ){ + + input=(Input*)(*object); + input->Configure(parameters); + + } + +} +/*}}}*/ Index: /issm/trunk/src/c/Container/Inputs.h =================================================================== --- /issm/trunk/src/c/Container/Inputs.h (revision 8362) +++ /issm/trunk/src/c/Container/Inputs.h (revision 8363) @@ -48,4 +48,6 @@ void GetParameterValue(int* pvalue,int enum_type); void GetParameterValue(double* pvalue,int enum_type); + + void Configure(Parameters* parameters); /*}}}*/ Index: /issm/trunk/src/c/EnumDefinitions/EnumDefinitions.h =================================================================== --- /issm/trunk/src/c/EnumDefinitions/EnumDefinitions.h (revision 8362) +++ /issm/trunk/src/c/EnumDefinitions/EnumDefinitions.h (revision 8363) @@ -135,4 +135,5 @@ PentaVertexInputEnum, TriaVertexInputEnum, + TriaVertexForcingEnum, ControlInputEnum, /*Params: */ Index: /issm/trunk/src/c/Makefile.am =================================================================== --- /issm/trunk/src/c/Makefile.am (revision 8362) +++ /issm/trunk/src/c/Makefile.am (revision 8363) @@ -194,4 +194,6 @@ ./objects/Inputs/TriaVertexInput.h\ ./objects/Inputs/TriaVertexInput.cpp\ + ./objects/Inputs/TriaVertexForcing.h\ + ./objects/Inputs/TriaVertexForcing.cpp\ ./objects/Inputs/PentaVertexInput.h\ ./objects/Inputs/PentaVertexInput.cpp\ @@ -839,4 +841,6 @@ ./objects/Inputs/TriaVertexInput.h\ ./objects/Inputs/TriaVertexInput.cpp\ + ./objects/Inputs/TriaVertexForcing.h\ + ./objects/Inputs/TriaVertexForcing.cpp\ ./objects/Inputs/PentaVertexInput.h\ ./objects/Inputs/PentaVertexInput.cpp\ Index: /issm/trunk/src/c/modules/EnumToStringx/EnumToStringx.cpp =================================================================== --- /issm/trunk/src/c/modules/EnumToStringx/EnumToStringx.cpp (revision 8362) +++ /issm/trunk/src/c/modules/EnumToStringx/EnumToStringx.cpp (revision 8363) @@ -114,4 +114,5 @@ case PentaVertexInputEnum : return "PentaVertexInput"; case TriaVertexInputEnum : return "TriaVertexInput"; + case TriaVertexForcingEnum : return "TriaVertexForcing"; case ControlInputEnum : return "ControlInput"; case ParamEnum : return "Param"; Index: /issm/trunk/src/c/modules/ModelProcessorx/CreateDataSets.cpp =================================================================== --- /issm/trunk/src/c/modules/ModelProcessorx/CreateDataSets.cpp (revision 8362) +++ /issm/trunk/src/c/modules/ModelProcessorx/CreateDataSets.cpp (revision 8363) @@ -115,5 +115,5 @@ /*Update Elements in case we are running a transient solution: */ - if(analysis_counter==nummodels && (solution_type==Transient2DSolutionEnum | solution_type==Transient3DSolutionEnum)){ + if(analysis_counter==(nummodels-1)&& (solution_type==Transient2DSolutionEnum || solution_type==Transient3DSolutionEnum)){ UpdateElementsTransient(elements,iomodel,iomodel_handle,analysis_counter,analysis_type); } Index: /issm/trunk/src/c/modules/ModelProcessorx/Transient/UpdateElementsTransient.cpp =================================================================== --- /issm/trunk/src/c/modules/ModelProcessorx/Transient/UpdateElementsTransient.cpp (revision 8362) +++ /issm/trunk/src/c/modules/ModelProcessorx/Transient/UpdateElementsTransient.cpp (revision 8363) @@ -16,10 +16,52 @@ /*Intermediary*/ - int i; + int i,j; int counter; Element *element = NULL; + char fetchstring[100]; + double time; + + + /*download whatever data will persist through forcing processing: */ + IoModelFetchData(&iomodel->elements,NULL,NULL,iomodel_handle,"elements"); /*Ok, let's go through forcings: */ + /*first, accumulation: */ + if(iomodel->forcing_accumulation_num_time_steps==0)_error_("accumulation forcing not available!"); + + iomodel->forcing_numtimesteps=iomodel->forcing_accumulation_num_time_steps; + IoModelFetchData(&iomodel->timesteps,NULL,NULL,iomodel_handle,"forcing_accumulation_time_steps"); + + for(i=0;iforcing_accumulation_num_time_steps;i++){ + time=iomodel->timesteps[i]; + sprintf(&fetchstring[0],"forcing_accumulation_%i",i+1); + IoModelFetchData(&iomodel->forcing,NULL,NULL,iomodel_handle,fetchstring); + + /*we now have the forcing for the corresponding time step, for all the nodes. Use this + *to write over the existing accumulation input: */ + + counter=0; + for (j=0;jnumberofelements;j++){ + if(iomodel->my_elements[j]){ + element=(Element*)elements->GetObjectByOffset(counter); + element->UpdateForcing(j,iomodel,i,time,AccumulationRateEnum); //we need j to index into elements. + counter++; + + } + } + + + + /*Free ressources:*/ + xfree((void**)&iomodel->forcing); + } + + + /*Free ressources:*/ + xfree((void**)&iomodel->forcing); + xfree((void**)&iomodel->timesteps); + xfree((void**)&iomodel->elements); + } Index: /issm/trunk/src/c/modules/StringToEnumx/StringToEnumx.cpp =================================================================== --- /issm/trunk/src/c/modules/StringToEnumx/StringToEnumx.cpp (revision 8362) +++ /issm/trunk/src/c/modules/StringToEnumx/StringToEnumx.cpp (revision 8363) @@ -112,4 +112,5 @@ else if (strcmp(name,"PentaVertexInput")==0) return PentaVertexInputEnum; else if (strcmp(name,"TriaVertexInput")==0) return TriaVertexInputEnum; + else if (strcmp(name,"TriaVertexForcing")==0) return TriaVertexForcingEnum; else if (strcmp(name,"ControlInput")==0) return ControlInputEnum; else if (strcmp(name,"Param")==0) return ParamEnum; Index: /issm/trunk/src/c/objects/Elements/Element.h =================================================================== --- /issm/trunk/src/c/objects/Elements/Element.h (revision 8362) +++ /issm/trunk/src/c/objects/Elements/Element.h (revision 8363) @@ -59,4 +59,5 @@ virtual void DeleteResults(void)=0; virtual void Update(int index, IoModel* iomodel,int analysis_counter,int analysis_type)=0; + virtual void UpdateForcing(int index, IoModel* iomodel,int step,double time, int FieldEnum)=0; virtual void InputToResult(int enum_type,int step,double time)=0; virtual void ControlInputGetGradient(Vec gradient,int enum_type)=0; Index: /issm/trunk/src/c/objects/Elements/Penta.h =================================================================== --- /issm/trunk/src/c/objects/Elements/Penta.h (revision 8362) +++ /issm/trunk/src/c/objects/Elements/Penta.h (revision 8363) @@ -129,4 +129,5 @@ double SurfaceArea(void); void Update(int index, IoModel* iomodel,int analysis_counter,int analysis_type); + void UpdateForcing(int index, IoModel* iomodel,int step,double time, int FieldEnum){_error_("not supported yet!");} int UpdateShelfStatus(Vec new_shelf_nodes); void UpdateShelfFlags(double* new_shelf_nodes); Index: /issm/trunk/src/c/objects/Elements/Tria.cpp =================================================================== --- /issm/trunk/src/c/objects/Elements/Tria.cpp (revision 8362) +++ /issm/trunk/src/c/objects/Elements/Tria.cpp (revision 8363) @@ -2725,4 +2725,8 @@ /*point parameters to real dataset: */ this->parameters=parametersin; + + /*get inputs configured too: */ + this->inputs->Configure(parameters); + } @@ -5863,4 +5867,37 @@ } /*}}}*/ +/*FUNCTION Tria::UpdateForcing(int index,IoModel* iomodel,int step,double time, int FieldEnum){{{1*/ +void Tria::UpdateForcing(int index,IoModel* iomodel,int step,double time, int FieldEnum){ + + /*Intermediaries*/ + int i,j; + int tria_vertex_ids[3]; + double nodeinputs[3]; + + /*Recover vertices ids needed to initialize inputs*/ + for(i=0;i<3;i++) tria_vertex_ids[i]=(int)iomodel->elements[3*index+i]; //ids for vertices are in the elements array from Matlab + + /*recover forcing for this time step: */ + for(i=0;i<3;i++)nodeinputs[i]=iomodel->forcing[tria_vertex_ids[i]-1]; + + if(step==0){ + /*This is the first time we are trying to replace the accumulation tria vertex input by an accumulation triav vertex + *forcing, which is essentially a tria vertex input that can vary with time. So firt, if there is already + an input with enum FieldEnum, squash it: */ + this->inputs->AddInput(new TriaVertexForcing(FieldEnum,nodeinputs,time,iomodel->forcing_numtimesteps,this->parameters)); + } + else{ + + /*Ok, we can't write over the existing forcing that was created in previous step 0. Find the input, and add + *values for this time step: */ + TriaVertexForcing* forcing=(TriaVertexForcing*)inputs->GetInput(FieldEnum); _assert_(forcing); + + forcing->AddTimeValues(&nodeinputs[0],step,time); + + } + + +} +/*}}}*/ /*FUNCTION Tria::UpdateShelfStatus{{{1*/ int Tria::UpdateShelfStatus(Vec new_shelf_nodes){ Index: /issm/trunk/src/c/objects/Elements/Tria.h =================================================================== --- /issm/trunk/src/c/objects/Elements/Tria.h (revision 8362) +++ /issm/trunk/src/c/objects/Elements/Tria.h (revision 8363) @@ -134,4 +134,5 @@ double SurfaceArea(void); void Update(int index, IoModel* iomodel,int analysis_counter,int analysis_type); + void UpdateForcing(int index, IoModel* iomodel,int step,double time, int FieldEnum); int UpdateShelfStatus(Vec new_shelf_nodes); int UpdatePotentialSheetUngrounding(double* potential_sheet_ungrounding,Vec vec_nodes_on_iceshelf,double* nodes_on_iceshelf); Index: /issm/trunk/src/c/objects/Inputs/BoolInput.cpp =================================================================== --- /issm/trunk/src/c/objects/Inputs/BoolInput.cpp (revision 8362) +++ /issm/trunk/src/c/objects/Inputs/BoolInput.cpp (revision 8363) @@ -237,2 +237,7 @@ } /*}}}*/ +/*FUNCTION BoolInput::Configure{{{1*/ +void BoolInput::Configure(Parameters* parameters){ + /*do nothing: */ +} +/*}}}*/ Index: /issm/trunk/src/c/objects/Inputs/BoolInput.h =================================================================== --- /issm/trunk/src/c/objects/Inputs/BoolInput.h (revision 8362) +++ /issm/trunk/src/c/objects/Inputs/BoolInput.h (revision 8363) @@ -44,4 +44,6 @@ Input* PointwiseMax(Input* inputB){_error_("not implemented yet");}; ElementResult* SpawnResult(int step, double time); + void Configure(Parameters* parameters); + void AddTimeValues(double* values,int step,double time){_error_("not supported yet");}; /*}}}*/ /*numerics: {{{1*/ Index: /issm/trunk/src/c/objects/Inputs/ControlInput.cpp =================================================================== --- /issm/trunk/src/c/objects/Inputs/ControlInput.cpp (revision 8362) +++ /issm/trunk/src/c/objects/Inputs/ControlInput.cpp (revision 8363) @@ -448,2 +448,7 @@ values->VerticallyIntegrate(thickness_input); }/*}}}*/ +/*FUNCTION ControlInput::Configure{{{1*/ +void ControlInput::Configure(Parameters* parameters){ + /*do nothing: */ +} +/*}}}*/ Index: /issm/trunk/src/c/objects/Inputs/ControlInput.h =================================================================== --- /issm/trunk/src/c/objects/Inputs/ControlInput.h (revision 8362) +++ /issm/trunk/src/c/objects/Inputs/ControlInput.h (revision 8363) @@ -48,4 +48,6 @@ Input* PointwiseMax(Input* inputB){_error_("not implemented yet");}; ElementResult* SpawnResult(int step, double time); + void AddTimeValues(double* values,int step,double time){_error_("not supported yet");}; + void Configure(Parameters* parameters); /*}}}*/ /*numerics: {{{1*/ Index: /issm/trunk/src/c/objects/Inputs/DoubleInput.cpp =================================================================== --- /issm/trunk/src/c/objects/Inputs/DoubleInput.cpp (revision 8362) +++ /issm/trunk/src/c/objects/Inputs/DoubleInput.cpp (revision 8363) @@ -368,2 +368,7 @@ } /*}}}*/ +/*FUNCTION DoubleInput::Configure{{{1*/ +void DoubleInput::Configure(Parameters* parameters){ + /*do nothing: */ +} +/*}}}*/ Index: /issm/trunk/src/c/objects/Inputs/DoubleInput.h =================================================================== --- /issm/trunk/src/c/objects/Inputs/DoubleInput.h (revision 8362) +++ /issm/trunk/src/c/objects/Inputs/DoubleInput.h (revision 8363) @@ -43,4 +43,6 @@ Input* PointwiseMax(Input* inputB); ElementResult* SpawnResult(int step, double time); + void AddTimeValues(double* values,int step,double time){_error_("not supported yet");}; + void Configure(Parameters* parameters); /*}}}*/ /*numerics: {{{1*/ Index: /issm/trunk/src/c/objects/Inputs/Input.h =================================================================== --- /issm/trunk/src/c/objects/Inputs/Input.h (revision 8362) +++ /issm/trunk/src/c/objects/Inputs/Input.h (revision 8363) @@ -13,4 +13,5 @@ class ElementResult; class GaussTria; +class Parameters; /*}}}*/ @@ -38,4 +39,5 @@ virtual void GetVyStrainRate3dPattyn(double* epsilonvy,double* xyz_list, GaussPenta* gauss)=0; virtual void ChangeEnum(int newenumtype)=0; + virtual void Configure(Parameters* parameters)=0; virtual void SquareMin(double* psquaremin, bool process_units,Parameters* parameters)=0; @@ -60,4 +62,5 @@ virtual Input* PointwiseMin(Input* inputmin)=0; virtual ElementResult* SpawnResult(int step, double time)=0; + virtual void AddTimeValues(double* values,int step,double time)=0; /*}}}*/ Index: /issm/trunk/src/c/objects/Inputs/IntInput.cpp =================================================================== --- /issm/trunk/src/c/objects/Inputs/IntInput.cpp (revision 8362) +++ /issm/trunk/src/c/objects/Inputs/IntInput.cpp (revision 8363) @@ -242,2 +242,7 @@ } /*}}}*/ +/*FUNCTION IntInput::Configure{{{1*/ +void IntInput::Configure(Parameters* parameters){ + /*do nothing: */ +} +/*}}}*/ Index: /issm/trunk/src/c/objects/Inputs/IntInput.h =================================================================== --- /issm/trunk/src/c/objects/Inputs/IntInput.h (revision 8362) +++ /issm/trunk/src/c/objects/Inputs/IntInput.h (revision 8363) @@ -44,4 +44,6 @@ Input* PointwiseMax(Input* inputB){_error_("not implemented yet");}; ElementResult* SpawnResult(int step, double time); + void AddTimeValues(double* values,int step,double time){_error_("not supported yet");}; + void Configure(Parameters* parameters); /*}}}*/ /*numerics: {{{1*/ Index: /issm/trunk/src/c/objects/Inputs/PentaVertexInput.cpp =================================================================== --- /issm/trunk/src/c/objects/Inputs/PentaVertexInput.cpp (revision 8362) +++ /issm/trunk/src/c/objects/Inputs/PentaVertexInput.cpp (revision 8363) @@ -673,2 +673,7 @@ } /*}}}*/ +/*FUNCTION PentaVertexInput::Configure{{{1*/ +void PentaVertexInput::Configure(Parameters* parameters){ + /*do nothing: */ +} +/*}}}*/ Index: /issm/trunk/src/c/objects/Inputs/PentaVertexInput.h =================================================================== --- /issm/trunk/src/c/objects/Inputs/PentaVertexInput.h (revision 8362) +++ /issm/trunk/src/c/objects/Inputs/PentaVertexInput.h (revision 8363) @@ -44,4 +44,6 @@ Input* PointwiseMax(Input* inputB); ElementResult* SpawnResult(int step, double time); + void AddTimeValues(double* values,int step,double time){_error_("not supported yet");}; + void Configure(Parameters* parameters); /*}}}*/ /*numerics: {{{1*/ Index: /issm/trunk/src/c/objects/Inputs/TriaVertexForcing.cpp =================================================================== --- /issm/trunk/src/c/objects/Inputs/TriaVertexForcing.cpp (revision 8363) +++ /issm/trunk/src/c/objects/Inputs/TriaVertexForcing.cpp (revision 8363) @@ -0,0 +1,507 @@ +/*!\file TriaVertexForcing.c + * \brief: implementation of the TriaVertexForcing object + */ + +#ifdef HAVE_CONFIG_H + #include "config.h" +#else +#error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!" +#endif + +#include "stdio.h" +#include +#include "../objects.h" +#include "../../EnumDefinitions/EnumDefinitions.h" +#include "../../shared/shared.h" +#include "../../Container/Container.h" +#include "../../include/include.h" + +/*TriaVertexForcing constructors and destructor*/ +/*FUNCTION TriaVertexForcing::TriaVertexForcing(){{{1*/ +TriaVertexForcing::TriaVertexForcing(){ + enum_type=NoneEnum; + values=NULL; + numtimesteps=0; + parameters=NULL; + return; +} +/*}}}*/ +/*FUNCTION TriaVertexForcing::TriaVertexForcing(int in_enum_type,double* time0values,double time0,int numtimesteps,Parameters* parameters){{{1*/ +TriaVertexForcing::TriaVertexForcing(int in_enum_type,double* time0values,double time0,int in_numtimesteps,Parameters* parameters) + :TriaRef(1) +{ + + /*Set TriaRef*/ + this->SetElementType(P1Enum,0); + this->element_type=P1Enum; + + /*Set Enum*/ + enum_type=in_enum_type; + + /*Allocate values and timesteps, and set first time step values: */ + this->numtimesteps=in_numtimesteps; + this->values=(double*)xmalloc(this->numtimesteps*3*sizeof(double)); + this->timesteps=(double*)xmalloc(this->numtimesteps*sizeof(double)); + this->values[0]=time0values[0]; + this->values[1]=time0values[1]; + this->values[2]=time0values[2]; + this->timesteps[0]=time0; + + this->parameters=parameters; + +} +/*}}}*/ +/*FUNCTION TriaVertexForcing::TriaVertexForcing(int in_enum_type,double* in_values,double* in_time,int numtimesteps,Parameters* parameters){{{1*/ +TriaVertexForcing::TriaVertexForcing(int in_enum_type,double* in_values,double* in_time,int in_numtimesteps,Parameters* parameters) + :TriaRef(1) +{ + + /*Set TriaRef*/ + this->SetElementType(P1Enum,0); + this->element_type=P1Enum; + + /*Set Enum*/ + enum_type=in_enum_type; + + /*Allocate values and timesteps, and copy: */ + this->numtimesteps=in_numtimesteps; + this->values=(double*)xmalloc(this->numtimesteps*3*sizeof(double)); + this->timesteps=(double*)xmalloc(this->numtimesteps*sizeof(double)); + + memcpy(this->values,in_values,numtimesteps*3*sizeof(double)); + memcpy(this->timesteps,in_time,numtimesteps*sizeof(double)); + + this->parameters=parameters; + +} +/*}}}*/ +/*FUNCTION TriaVertexForcing::~TriaVertexForcing(){{{1*/ +TriaVertexForcing::~TriaVertexForcing(){ + xfree((void**)&this->values); + xfree((void**)&this->timesteps); + parameters=NULL; + return; +} +/*}}}*/ +/*FUNCTION TriaVertexForcing::AddTimeValues(double* values,int step,double time);{{{1*/ +void TriaVertexForcing::AddTimeValues(double* values,int step,double time){ + + /*insert values at time step: */ + if (step<0)_error_("timestep should start at least at time 0!"); + if (step>numtimesteps-1)_error_("timestep cannot be more than number of time steps"); + + /*go ahead and plug: */ + this->timesteps[step]=time; + this->values[3*step+0]=values[0]; + this->values[3*step+1]=values[1]; + this->values[3*step+2]=values[2]; +} +/*}}}*/ + +/*Object virtual functions definitions:*/ +/*FUNCTION TriaVertexForcing::Echo {{{1*/ +void TriaVertexForcing::Echo(void){ + this->DeepEcho(); +} +/*}}}*/ +/*FUNCTION TriaVertexForcing::DeepEcho{{{1*/ +void TriaVertexForcing::DeepEcho(void){ + + int i; + + printf("TriaVertexForcing:\n"); + printf(" enum: %i (%s)\n",this->enum_type,EnumToStringx(this->enum_type)); + printf(" numtimesteps: %i\n",this->numtimesteps); + for(i=0;inumtimesteps;i++){ + printf(" time: %g values: [%g %g %g]\n",this->timesteps[i],this->values[3*i+0],this->values[3*i+1],this->values[3*i+2]); + } +} +/*}}}*/ +/*FUNCTION TriaVertexForcing::Id{{{1*/ +int TriaVertexForcing::Id(void){ return -1; } +/*}}}*/ +/*FUNCTION TriaVertexForcing::MyRank{{{1*/ +int TriaVertexForcing::MyRank(void){ + extern int my_rank; + return my_rank; +} +/*}}}*/ +/*FUNCTION TriaVertexForcing::Marshall{{{1*/ +void TriaVertexForcing::Marshall(char** pmarshalled_dataset){ + + char* marshalled_dataset=NULL; + int enum_value=0; + + /*recover marshalled_dataset: */ + marshalled_dataset=*pmarshalled_dataset; + + /*get enum value of TriaVertexForcing: */ + enum_value=TriaVertexForcingEnum; + + /*marshall enum: */ + memcpy(marshalled_dataset,&enum_value,sizeof(enum_value));marshalled_dataset+=sizeof(enum_value); + + /*marshall TriaVertexForcing data: */ + memcpy(marshalled_dataset,&enum_type,sizeof(enum_type));marshalled_dataset+=sizeof(enum_type); + memcpy(marshalled_dataset,&numtimesteps,sizeof(numtimesteps));marshalled_dataset+=sizeof(numtimesteps); + memcpy(marshalled_dataset,values,numtimesteps*3*sizeof(double));marshalled_dataset+=numtimesteps*3*sizeof(double); + memcpy(marshalled_dataset,timesteps,numtimesteps*sizeof(double));marshalled_dataset+=numtimesteps*sizeof(double); + + *pmarshalled_dataset=marshalled_dataset; +} +/*}}}*/ +/*FUNCTION TriaVertexForcing::MarshallSize{{{1*/ +int TriaVertexForcing::MarshallSize(){ + + return + +sizeof(enum_type)+ + +sizeof(numtimesteps)+ + 3*numtimesteps*sizeof(double)+ + numtimesteps*sizeof(double)+ + +sizeof(int); //sizeof(int) for enum value +} +/*}}}*/ +/*FUNCTION TriaVertexForcing::Demarshall{{{1*/ +void TriaVertexForcing::Demarshall(char** pmarshalled_dataset){ + + char* marshalled_dataset=NULL; + int i; + + /*recover marshalled_dataset: */ + marshalled_dataset=*pmarshalled_dataset; + + /*this time, no need to get enum type, the pointer directly points to the beginning of the + *object data (thanks to DataSet::Demarshall):*/ + memcpy(&enum_type,marshalled_dataset,sizeof(enum_type));marshalled_dataset+=sizeof(enum_type); + memcpy(&numtimesteps,marshalled_dataset,sizeof(numtimesteps));marshalled_dataset+=sizeof(numtimesteps); + /*allocate: */ + timesteps=(double*)xmalloc(numtimesteps*sizeof(double)); + values=(double*)xmalloc(numtimesteps*3*sizeof(double)); + + memcpy(values,marshalled_dataset,3*numtimesteps*sizeof(double));marshalled_dataset+=3*numtimesteps*sizeof(double); + memcpy(timesteps,marshalled_dataset,3*numtimesteps*sizeof(double));marshalled_dataset+=3*numtimesteps*sizeof(double); + + /*return: */ + *pmarshalled_dataset=marshalled_dataset; + return; +} +/*}}}*/ +/*FUNCTION TriaVertexForcing::Enum{{{1*/ +int TriaVertexForcing::Enum(void){ + + return TriaVertexForcingEnum; + +} +/*}}}*/ +/*FUNCTION TriaVertexForcing::copy{{{1*/ +Object* TriaVertexForcing::copy() { + + return new TriaVertexForcing(this->enum_type,this->values,this->timesteps,this->numtimesteps,this->parameters); + +} +/*}}}*/ + +/*TriaVertexForcing management*/ +/*FUNCTION TriaVertexForcing::EnumType{{{1*/ +int TriaVertexForcing::EnumType(void){ + + return this->enum_type; + +} +/*}}}*/ +/*FUNCTION TriaVertexForcing::SpawnTriaInput{{{1*/ +Input* TriaVertexForcing::SpawnTriaInput(int* indices){ + + /*output*/ + TriaVertexForcing* outforcing=NULL; + + /*Create new Tria forcing (copy of current forcing)*/ + outforcing=new TriaVertexForcing(this->enum_type,this->values,this->timesteps,this->numtimesteps,this->parameters); + + /*Assign output*/ + return outforcing; + +} +/*}}}*/ +/*FUNCTION TriaVertexForcing::SpawnResult{{{1*/ +ElementResult* TriaVertexForcing::SpawnResult(int step, double time){ + + double triavalues[3]; + + /*Ok, we want to spawn a TriaVertexElementResult. We have the time, just get + *the correct values at the three nodes: */ + this->GetTimeValues(&triavalues[0],time); + + return new TriaVertexElementResult(this->enum_type,&triavalues[0],step,time); + +} +/*}}}*/ + + +/*Object functions*/ +/*FUNCTION TriaVertexForcing::GetParameterValue(double* pvalue,GaussTria* gauss){{{1*/ +void TriaVertexForcing::GetParameterValue(double* pvalue,GaussTria* gauss){ + + double time; + double triavalues[3]; + + /*First, recover current time from parameters: */ + parameters->FindParam(&time,TimeEnum); + + /*Retrieve interpolated values for this time step: */ + this->GetTimeValues(&triavalues[0],time); + + /*Call TriaRef function*/ + TriaRef::GetParameterValue(pvalue,&triavalues[0],gauss); + +} +/*}}}*/ +/*FUNCTION TriaVertexForcing::GetParameterDerivativeValue(double* p, double* xyz_list, GaussTria* gauss){{{1*/ +void TriaVertexForcing::GetParameterDerivativeValue(double* p, double* xyz_list, GaussTria* gauss){ + + double time; + double triavalues[3]; + + /*First, recover current time from parameters: */ + parameters->FindParam(&time,TimeEnum); + + /*Retrieve interpolated values for this time step: */ + this->GetTimeValues(&triavalues[0],time); + + /*Call TriaRef function*/ + TriaRef::GetParameterDerivativeValue(p,&triavalues[0],xyz_list,gauss); +} +/*}}}*/ +/*FUNCTION TriaVertexForcing::ChangeEnum{{{1*/ +void TriaVertexForcing::ChangeEnum(int newenumtype){ + this->enum_type=newenumtype; +} +/*}}}*/ +/*FUNCTION TriaVertexForcing::GetParameterAverage{{{1*/ +void TriaVertexForcing::GetParameterAverage(double* pvalue){ + + double time; + double triavalues[3]; + + /*First, recover current time from parameters: */ + parameters->FindParam(&time,TimeEnum); + + /*Retrieve interpolated values for this time step: */ + this->GetTimeValues(&triavalues[0],time); + + *pvalue=1./3.*(triavalues[0]+triavalues[1]+triavalues[2]); +} +/*}}}*/ + +/*Intermediary*/ +/*FUNCTION TriaVertexForcing::SquareMin{{{1*/ +void TriaVertexForcing::SquareMin(double* psquaremin, bool process_units,Parameters* parameters){ + + int i; + const int numnodes=3; + double valuescopy[numnodes]; + double squaremin; + double time; + + /*First, recover current time from parameters: */ + parameters->FindParam(&time,TimeEnum); + + /*Retrieve interpolated values for this time step: */ + this->GetTimeValues(&valuescopy[0],time); + + /*Process units if requested: */ + if(process_units)UnitConversion(&valuescopy[0],numnodes,IuToExtEnum,enum_type,parameters); + + /*Now, figure out minimum of valuescopy: */ + squaremin=pow(valuescopy[0],2); + for(i=1;iFindParam(&time,TimeEnum); + + /*Retrieve interpolated values for this time step: */ + this->GetTimeValues(&triavalues[0],time); + + for(int i=0;inorm) norm=fabs(triavalues[i]); + return norm; +} +/*}}}*/ +/*FUNCTION TriaVertexForcing::Max{{{1*/ +double TriaVertexForcing::Max(void){ + + const int numnodes=3; + double max; + double time; + double triavalues[3]; + + /*First, recover current time from parameters: */ + parameters->FindParam(&time,TimeEnum); + + /*Retrieve interpolated values for this time step: */ + this->GetTimeValues(&triavalues[0],time); + + max=triavalues[0]; + + for(int i=1;imax) max=triavalues[i]; + } + return max; +} +/*}}}*/ +/*FUNCTION TriaVertexForcing::MaxAbs{{{1*/ +double TriaVertexForcing::MaxAbs(void){ + + const int numnodes=3; + double max; + double time; + double triavalues[3]; + + /*First, recover current time from parameters: */ + parameters->FindParam(&time,TimeEnum); + + /*Retrieve interpolated values for this time step: */ + this->GetTimeValues(&triavalues[0],time); + + + max=fabs(triavalues[0]); + + for(int i=1;imax) max=fabs(triavalues[i]); + } + return max; +} +/*}}}*/ +/*FUNCTION TriaVertexForcing::Min{{{1*/ +double TriaVertexForcing::Min(void){ + + const int numnodes=3; + double min; + double time; + double triavalues[3]; + + /*First, recover current time from parameters: */ + parameters->FindParam(&time,TimeEnum); + + /*Retrieve interpolated values for this time step: */ + this->GetTimeValues(&triavalues[0],time); + + min=triavalues[0]; + + for(int i=1;iFindParam(&time,TimeEnum); + + /*Retrieve interpolated values for this time step: */ + this->GetTimeValues(&triavalues[0],time); + + min=fabs(triavalues[0]); + + for(int i=1;iFindParam(&time,TimeEnum); + + /*Retrieve interpolated values for this time step: */ + this->GetTimeValues(&triavalues[0],time); + + VecSetValues(vector,numvertices,doflist,(const double*)&triavalues[0],INSERT_VALUES); + +} /*}}}*/ +/*FUNCTION TriaVertexForcing::GetTimeValues{{{1*/ +void TriaVertexForcing::GetTimeValues(double* outvalues,double intime){ + + int i,j; + double deltat; + double alpha; + bool found=false; + + /*Ok, we have the time, go through the timesteps, and figure out which interval we + *fall within. Then interpolate the values on this interval: */ + if(intimetimesteps[0]){ + /*get values for the first time: */ + outvalues[0]=this->values[3*0+0]; + outvalues[1]=this->values[3*0+1]; + outvalues[2]=this->values[3*0+2]; + found=true; + } + else if(intime>this->timesteps[this->numtimesteps-1]){ + /*get values for the last time: */ + + outvalues[0]=this->values[3*this->numtimesteps-1+0]; + outvalues[1]=this->values[3*this->numtimesteps-1+1]; + outvalues[2]=this->values[3*this->numtimesteps-1+2]; + found=true; + } + else{ + /*Find which interval we fall within: */ + for(i=0;inumtimesteps;i++){ + if(intime==this->timesteps[i]){ + /*We are right on one step time: */ + outvalues[0]=this->values[3*i+0]; + outvalues[1]=this->values[3*i+1]; + outvalues[2]=this->values[3*i+2]; + found=true; + break; //we are done with the time interpolation. + } + else{ + if(this->timesteps[i]timesteps[i+1]){ + /*ok, we have the interval ]i:i+1[. Interpolate linearly for now: */ + deltat=this->timesteps[i+1]-this->timesteps[i]; + alpha=(intime-this->timesteps[i])/deltat; + for(j=0;j<3;j++){ + outvalues[j]=this->values[3*i+j]+alpha*(this->values[3*(i+1)+j]-this->values[3*i+j]); + } + found=true; + break; + } + else continue; //keep looking on the next interval + } + } + } + if(!found)_error_("did not find time interval on which to interpolate forcing values!"); + +} +/*}}}*/ +/*FUNCTION TriaVertexForcing::Configure{{{1*/ +void TriaVertexForcing::Configure(Parameters* parameters){ + this->parameters=parameters; +} +/*}}}*/ Index: /issm/trunk/src/c/objects/Inputs/TriaVertexForcing.h =================================================================== --- /issm/trunk/src/c/objects/Inputs/TriaVertexForcing.h (revision 8363) +++ /issm/trunk/src/c/objects/Inputs/TriaVertexForcing.h (revision 8363) @@ -0,0 +1,91 @@ +/*! \file TriaVertexForcing.h + * \brief: header file for triavertexforcing object + */ + + +#ifndef _TRIAVERTEXFORCING_H_ +#define _TRIAVERTEXFORCING_H_ + +/*Headers:*/ +/*{{{1*/ +#include "./Input.h" +#include "../Elements/TriaRef.h" +class GaussTria; +class Parameters; +/*}}}*/ + +class TriaVertexForcing: public Input,public TriaRef{ + + public: + /*just hold 3 values for 3 vertices, for all time steps: */ + int enum_type; + int numtimesteps; + double* values; + double* timesteps; + Parameters* parameters; //to find current time. + + /*TriaVertexForcing constructors, destructors: {{{1*/ + TriaVertexForcing(); + TriaVertexForcing(int enum_type,double* valuestime0,double time0,int numtimesteps,Parameters* parameters); + TriaVertexForcing(int enum_type,double* values,double* timesteps,int numtimesteps,Parameters* parameters); + ~TriaVertexForcing(); + void AddTimeValues(double* values,int step,double time); + /*}}}*/ + /*Object virtual functions definitions:{{{1 */ + void Echo(); + void DeepEcho(); + int Id(); + int MyRank(); + void Marshall(char** pmarshalled_dataset); + int MarshallSize(); + void Demarshall(char** pmarshalled_dataset); + int Enum(); + Object* copy(); + /*}}}*/ + /*TriaVertexForcing management: {{{1*/ + int EnumType(); + Input* SpawnTriaInput(int* indices); + Input* PointwiseDivide(Input* forcingB){_error_("not implemented yet");}; + Input* PointwiseMin(Input* forcingB){_error_("not implemented yet");}; + Input* PointwiseMax(Input* forcingB){_error_("not implemented yet");}; + ElementResult* SpawnResult(int step, double time); + void Configure(Parameters* parameters); + /*}}}*/ + /*numerics: {{{1*/ + void GetParameterValue(bool* pvalue){_error_("not implemented yet");} + void GetParameterValue(int* pvalue){_error_("not implemented yet");} + void GetParameterValue(double* pvalue){_error_("not implemented yet");} + void GetParameterValue(double* pvalue,GaussTria* gauss); + void GetParameterValue(double* pvalue,GaussPenta* gauss){_error_("not implemented yet");}; + void GetParameterDerivativeValue(double* derivativevalues, double* xyz_list, GaussTria* gauss); + void GetParameterDerivativeValue(double* derivativevalues, double* xyz_list, GaussPenta* gauss){_error_("not implemented yet");}; + void GetParameterAverage(double* pvalue); + void GetVxStrainRate2d(double* epsilonvx,double* xyz_list, GaussTria* gauss){_error_("not implemented yet");}; + void GetVyStrainRate2d(double* epsilonvy,double* xyz_list, GaussTria* gauss){_error_("not implemented yet");}; + void GetVxStrainRate3d(double* epsilonvx,double* xyz_list, GaussPenta* gauss){_error_("not implemented yet");}; + void GetVyStrainRate3d(double* epsilonvy,double* xyz_list, GaussPenta* gauss){_error_("not implemented yet");}; + void GetVzStrainRate3d(double* epsilonvz,double* xyz_list, GaussPenta* gauss){_error_("not implemented yet");}; + void GetVxStrainRate3dPattyn(double* epsilonvx,double* xyz_list, GaussPenta* gauss){_error_("not implemented yet");}; + void GetVyStrainRate3dPattyn(double* epsilonvy,double* xyz_list, GaussPenta* gauss){_error_("not implemented yet");}; + void ChangeEnum(int newenumtype); + + void SquareMin(double* psquaremin, bool process_units,Parameters* parameters); + void ConstrainMin(double minimum){_error_("not implemented yet");}; + void Scale(double scale_factor){_error_("not implemented yet");}; + void ArtificialNoise(double min,double max){_error_("not implemented yet");}; + void AXPY(Input* xforcing,double scalar){_error_("not implemented yet");}; + void Constrain(double cm_min, double cm_max){_error_("not implemented yet");}; + double InfinityNorm(void); + double Max(void); + double MaxAbs(void); + double Min(void); + double MinAbs(void); + void Extrude(void){_error_("not supported yet");}; + void VerticallyIntegrate(Input* thickness_forcing){_error_("not supported yet");}; + void GetVectorFromInputs(Vec vector,int* doflist); + void GetValuesPtr(double** pvalues,int* pnum_values){_error_("not supported yet");}; + void GetTimeValues(double* values,double time); + /*}}}*/ + +}; +#endif /* _TRIAVERTEXFORCING_H */ Index: /issm/trunk/src/c/objects/Inputs/TriaVertexInput.cpp =================================================================== --- /issm/trunk/src/c/objects/Inputs/TriaVertexInput.cpp (revision 8362) +++ /issm/trunk/src/c/objects/Inputs/TriaVertexInput.cpp (revision 8363) @@ -465,2 +465,7 @@ } /*}}}*/ +/*FUNCTION TriaVertexInput::Configure{{{1*/ +void TriaVertexInput::Configure(Parameters* parameters){ + /*do nothing: */ +} +/*}}}*/ Index: /issm/trunk/src/c/objects/Inputs/TriaVertexInput.h =================================================================== --- /issm/trunk/src/c/objects/Inputs/TriaVertexInput.h (revision 8362) +++ /issm/trunk/src/c/objects/Inputs/TriaVertexInput.h (revision 8363) @@ -44,4 +44,6 @@ Input* PointwiseMax(Input* inputB); ElementResult* SpawnResult(int step, double time); + void AddTimeValues(double* values,int step,double time){_error_("not supported yet");}; + void Configure(Parameters* parameters); /*}}}*/ /*numerics: {{{1*/ Index: /issm/trunk/src/c/objects/IoModel.cpp =================================================================== --- /issm/trunk/src/c/objects/IoModel.cpp (revision 8362) +++ /issm/trunk/src/c/objects/IoModel.cpp (revision 8363) @@ -89,5 +89,6 @@ xfree((void**)&this->melting_rate_correction); xfree((void**)&this->accumulation_rate); - xfree((void**)&this->forcing_accumulation); + xfree((void**)&this->forcing); + xfree((void**)&this->timesteps); xfree((void**)&this->dhdt); xfree((void**)&this->rheology_B); @@ -220,4 +221,5 @@ IoModelFetchData(&this->gl_melting_rate,iomodel_handle,"gl_melting_rate"); IoModelFetchData(&this->rheology_law,iomodel_handle,"rheology_law"); + IoModelFetchData(&this->forcing_accumulation_num_time_steps,iomodel_handle,"forcing_accumulation_num_time_steps"); /*qmu: */ @@ -397,6 +399,11 @@ /*!basal: */ this->accumulation_rate=NULL; - this->forcing_accumulation=NULL; this->dhdt=NULL; + + /*forcings: */ + this->forcing_accumulation_num_time_steps=0; + this->forcing_numtimesteps=0; + this->forcing=NULL; + this->timesteps=NULL; /*parameter output: */ Index: /issm/trunk/src/c/objects/IoModel.h =================================================================== --- /issm/trunk/src/c/objects/IoModel.h (revision 8362) +++ /issm/trunk/src/c/objects/IoModel.h (revision 8363) @@ -195,6 +195,11 @@ int melting_rate_correction_apply; double* accumulation_rate; - double* forcing_accumulation; double* dhdt; + + /*forcings: */ + int forcing_accumulation_num_time_steps; + int forcing_numtimesteps; + double* forcing; + double* timesteps; /*parameter output: */ Index: /issm/trunk/src/c/objects/objects.h =================================================================== --- /issm/trunk/src/c/objects/objects.h (revision 8362) +++ /issm/trunk/src/c/objects/objects.h (revision 8363) @@ -80,4 +80,5 @@ #include "./Inputs/PentaVertexInput.h" #include "./Inputs/TriaVertexInput.h" +#include "./Inputs/TriaVertexForcing.h" #include "./Inputs/ControlInput.h" Index: /issm/trunk/src/c/solutions/transient2d_core.cpp =================================================================== --- /issm/trunk/src/c/solutions/transient2d_core.cpp (revision 8362) +++ /issm/trunk/src/c/solutions/transient2d_core.cpp (revision 8363) @@ -56,4 +56,5 @@ } time+=dt; + femmodel->parameters->SetParam(time,TimeEnum); step+=1; @@ -85,4 +86,5 @@ InputToResultx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,SurfaceEnum,step,time); InputToResultx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,BedEnum,step,time); + InputToResultx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,AccumulationRateEnum,step,time); if(gl_migration!=NoneEnum)InputToResultx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,ElementOnIceShelfEnum,step,time); Index: /issm/trunk/src/c/solutions/transient3d_core.cpp =================================================================== --- /issm/trunk/src/c/solutions/transient3d_core.cpp (revision 8362) +++ /issm/trunk/src/c/solutions/transient3d_core.cpp (revision 8363) @@ -53,4 +53,5 @@ step+=1; time+=dt; + femmodel->parameters->SetParam(time,TimeEnum); _printf_(VerboseSolution(),"%s%g%s%i%s%g%s%g\n","time [yr]: ",time/yts," iteration number: ",step,"/",floor((finaltime-time)/dt)," dt [yr]: ",dt/yts); Index: /issm/trunk/src/m/enum/TriaVertexForcingEnum.m =================================================================== --- /issm/trunk/src/m/enum/TriaVertexForcingEnum.m (revision 8363) +++ /issm/trunk/src/m/enum/TriaVertexForcingEnum.m (revision 8363) @@ -0,0 +1,11 @@ +function macro=TriaVertexForcingEnum() +%TRIAVERTEXFORCINGENUM - Enum of TriaVertexForcing +% +% WARNING: DO NOT MODIFY THIS FILE +% this file has been automatically generated by src/c/EnumDefinitions/Synchronize.sh +% Please read src/c/EnumDefinitions/README for more information +% +% Usage: +% macro=TriaVertexForcingEnum() + +macro=StringToEnum('TriaVertexForcing'); Index: /issm/trunk/src/m/model/marshall.m =================================================================== --- /issm/trunk/src/m/model/marshall.m (revision 8362) +++ /issm/trunk/src/m/model/marshall.m (revision 8363) @@ -93,5 +93,15 @@ WriteData(fid,md.melting_rate,'Mat','melting_rate'); WriteData(fid,md.melting_rate_correction_apply,'Integer','melting_rate_correction_apply'); -WriteData(fid,md.forcing_accumulation,'Mat','forcing_accumulation'); + +if ~isnan(md.forcing_accumulation), + WriteData(fid,size(md.forcing_accumulation,2),'Integer','forcing_accumulation_num_time_steps'); + WriteData(fid,md.forcing_accumulation(end,:),'Mat','forcing_accumulation_time_steps'); + for i=1:size(md.forcing_accumulation,2), + WriteData(fid,md.forcing_accumulation(1:end-1,i),'Mat',['forcing_accumulation_' num2str(i)]); + end +else + WriteData(fid,0,'Integer','forcing_accumulation_num_time_steps'); +end + if md.melting_rate_correction_apply, WriteData(fid,md.melting_rate_correction,'Mat','melting_rate_correction'); Index: /issm/trunk/src/m/solutions/transient2d_core.m =================================================================== --- /issm/trunk/src/m/solutions/transient2d_core.m (revision 8362) +++ /issm/trunk/src/m/solutions/transient2d_core.m (revision 8363) @@ -33,4 +33,5 @@ step=step+1; time=time+dt; + femmodel.parameters.Time=time; issmprintf(VerboseSolution,'\n%s%g%s%i%s%g\n','time [yr]: ',time,' iteration number: ',step,'/',floor(ndt/dt)); Index: /issm/trunk/src/m/solutions/transient3d_core.m =================================================================== --- /issm/trunk/src/m/solutions/transient3d_core.m (revision 8362) +++ /issm/trunk/src/m/solutions/transient3d_core.m (revision 8363) @@ -32,4 +32,5 @@ step=step+1; time=time+dt; + femmodel.parameters.Time=time; issmprintf(VerboseSolution,'\n%s%g%s%i%s%g\n','time [yr] ',time,' iteration number: ',step,'/',floor(ndt/dt));