Index: /issm/trunk/src/c/Container/Parameters.cpp =================================================================== --- /issm/trunk/src/c/Container/Parameters.cpp (revision 8375) +++ /issm/trunk/src/c/Container/Parameters.cpp (revision 8376) @@ -57,5 +57,5 @@ /*}}}*/ /*FUNCTION Parameters::FindParam(bool* pbool,int enum_type){{{1*/ -void Parameters::FindParam(bool* pbool,int enum_type){ +void Parameters::FindParam(bool* pbool,int enum_type){ _assert_(this); vector::iterator object; @@ -74,5 +74,5 @@ /*}}}*/ /*FUNCTION Parameters::FindParam(int* pinteger,int enum_type){{{1*/ -void Parameters::FindParam(int* pinteger,int enum_type){ +void Parameters::FindParam(int* pinteger,int enum_type){ _assert_(this); vector::iterator object; @@ -91,5 +91,5 @@ /*}}}*/ /*FUNCTION Parameters::FindParam(double* pscalar, int enum_type){{{1*/ -void Parameters::FindParam(double* pscalar, int enum_type){ +void Parameters::FindParam(double* pscalar, int enum_type){ _assert_(this); vector::iterator object; @@ -108,5 +108,5 @@ /*}}}*/ /*FUNCTION Parameters::FindParam(char** pstring,int enum_type){{{1*/ -void Parameters::FindParam(char** pstring,int enum_type){ +void Parameters::FindParam(char** pstring,int enum_type){ _assert_(this); vector::iterator object; @@ -126,5 +126,5 @@ /*}}}*/ /*FUNCTION Parameters::FindParam(char*** pstringarray,int* pM,int enum_type){{{1*/ -void Parameters::FindParam(char*** pstringarray,int* pM,int enum_type){ +void Parameters::FindParam(char*** pstringarray,int* pM,int enum_type){ _assert_(this); vector::iterator object; @@ -144,5 +144,5 @@ /*}}}*/ /*FUNCTION Parameters::FindParam(int** pintarray,int* pM,int enum_type){{{1*/ -void Parameters::FindParam(int** pintarray,int* pM, int enum_type){ +void Parameters::FindParam(int** pintarray,int* pM, int enum_type){ _assert_(this); vector::iterator object; @@ -162,5 +162,5 @@ /*}}}*/ /*FUNCTION Parameters::FindParam(double** pdoublearray,int* pM,int enum_type){{{1*/ -void Parameters::FindParam(double** pdoublearray,int* pM, int enum_type){ +void Parameters::FindParam(double** pdoublearray,int* pM, int enum_type){ _assert_(this); vector::iterator object; @@ -180,5 +180,5 @@ /*}}}*/ /*FUNCTION Parameters::FindParam(double** pdoublearray,int* pM, int* pN,int enum_type){{{1*/ -void Parameters::FindParam(double** pdoublearray,int* pM, int* pN,int enum_type){ +void Parameters::FindParam(double** pdoublearray,int* pM, int* pN,int enum_type){ _assert_(this); vector::iterator object; @@ -198,5 +198,5 @@ /*}}}*/ /*FUNCTION Parameters::FindParam(double*** parray,int* pM,int** pmdims_array,int** pndims_array,int enum_type){{{1*/ -void Parameters::FindParam(double*** parray,int* pM,int** pmdims_array,int** pndims_array,int enum_type){ +void Parameters::FindParam(double*** parray,int* pM,int** pmdims_array,int** pndims_array,int enum_type){ _assert_(this); vector::iterator object; @@ -215,5 +215,5 @@ /*}}}*/ /*FUNCTION Parameters::FindParam(Vec* pvec,int enum_type){{{1*/ -void Parameters::FindParam(Vec* pvec,int enum_type){ +void Parameters::FindParam(Vec* pvec,int enum_type){ _assert_(this); vector::iterator object; @@ -233,5 +233,5 @@ /*}}}*/ /*FUNCTION Parameters::FindParam(Mat* pmat,int enum_type){{{1*/ -void Parameters::FindParam(Mat* pmat,int enum_type){ +void Parameters::FindParam(Mat* pmat,int enum_type){ _assert_(this); vector::iterator object; @@ -251,5 +251,5 @@ /*}}}*/ /*FUNCTION Parameters::FindParam(FILE** pfid,int enum_type){{{1*/ -void Parameters::FindParam(FILE** pfid,int enum_type){ +void Parameters::FindParam(FILE** pfid,int enum_type){ _assert_(this); vector::iterator object; Index: /issm/trunk/src/c/EnumDefinitions/EnumDefinitions.h =================================================================== --- /issm/trunk/src/c/EnumDefinitions/EnumDefinitions.h (revision 8375) +++ /issm/trunk/src/c/EnumDefinitions/EnumDefinitions.h (revision 8376) @@ -136,4 +136,5 @@ TriaVertexInputEnum, TriaVertexForcingEnum, + PentaVertexForcingEnum, ControlInputEnum, /*Params: */ Index: /issm/trunk/src/c/Makefile.am =================================================================== --- /issm/trunk/src/c/Makefile.am (revision 8375) +++ /issm/trunk/src/c/Makefile.am (revision 8376) @@ -198,4 +198,6 @@ ./objects/Inputs/PentaVertexInput.h\ ./objects/Inputs/PentaVertexInput.cpp\ + ./objects/Inputs/PentaVertexForcing.h\ + ./objects/Inputs/PentaVertexForcing.cpp\ ./objects/Inputs/BoolInput.h\ ./objects/Inputs/BoolInput.cpp\ @@ -845,4 +847,6 @@ ./objects/Inputs/PentaVertexInput.h\ ./objects/Inputs/PentaVertexInput.cpp\ + ./objects/Inputs/PentaVertexForcing.h\ + ./objects/Inputs/PentaVertexForcing.cpp\ ./objects/Inputs/BoolInput.h\ ./objects/Inputs/BoolInput.cpp\ Index: /issm/trunk/src/c/modules/EnumToStringx/EnumToStringx.cpp =================================================================== --- /issm/trunk/src/c/modules/EnumToStringx/EnumToStringx.cpp (revision 8375) +++ /issm/trunk/src/c/modules/EnumToStringx/EnumToStringx.cpp (revision 8376) @@ -115,4 +115,5 @@ case TriaVertexInputEnum : return "TriaVertexInput"; case TriaVertexForcingEnum : return "TriaVertexForcing"; + case PentaVertexForcingEnum : return "PentaVertexForcing"; case ControlInputEnum : return "ControlInput"; case ParamEnum : return "Param"; Index: /issm/trunk/src/c/modules/StringToEnumx/StringToEnumx.cpp =================================================================== --- /issm/trunk/src/c/modules/StringToEnumx/StringToEnumx.cpp (revision 8375) +++ /issm/trunk/src/c/modules/StringToEnumx/StringToEnumx.cpp (revision 8376) @@ -113,4 +113,5 @@ else if (strcmp(name,"TriaVertexInput")==0) return TriaVertexInputEnum; else if (strcmp(name,"TriaVertexForcing")==0) return TriaVertexForcingEnum; + else if (strcmp(name,"PentaVertexForcing")==0) return PentaVertexForcingEnum; else if (strcmp(name,"ControlInput")==0) return ControlInputEnum; else if (strcmp(name,"Param")==0) return ParamEnum; Index: /issm/trunk/src/c/objects/Elements/Penta.cpp =================================================================== --- /issm/trunk/src/c/objects/Elements/Penta.cpp (revision 8375) +++ /issm/trunk/src/c/objects/Elements/Penta.cpp (revision 8376) @@ -467,4 +467,7 @@ /*point parameters to real dataset: */ this->parameters=parametersin; + + /*get inputs configured too: */ + this->inputs->Configure(parameters); } /*}}}*/ @@ -6263,4 +6266,39 @@ } /*}}}*/ +/*FUNCTION Penta::UpdateForcing(int index,IoModel* iomodel,int step,double time, int FieldEnum){{{1*/ +void Penta::UpdateForcing(int index,IoModel* iomodel,int step,double time, int FieldEnum,Parameters* parameters){ + + /*Intermediaries*/ + int i,j; + int penta_vertex_ids[6]; + double nodeinputs[6]; + + /*Recover vertices ids needed to initialize inputs*/ + for(i=0;i<6;i++) penta_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<6;i++)nodeinputs[i]=iomodel->forcing[penta_vertex_ids[i]-1]; + + /*process units: */ + UnitConversion(&nodeinputs[0],6,ExtToIuEnum, FieldEnum, parameters); + + if(step==0){ + /*This is the first time we are trying to replace the accumulation penta vertex input by an accumulation pentav vertex + *forcing, which is essentially a penta vertex input that can vary with time. So firt, if there is already + an input with enum FieldEnum, squash it: */ + this->inputs->AddInput(new PentaVertexForcing(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: */ + PentaVertexForcing* forcing=(PentaVertexForcing*)inputs->GetInput(FieldEnum); _assert_(forcing); + forcing->AddTimeValues(&nodeinputs[0],step,time); + + } + + +} +/*}}}*/ /*FUNCTION Penta::UpdateShelfStatus{{{1*/ int Penta::UpdateShelfStatus(Vec new_shelf_nodes){ Index: /issm/trunk/src/c/objects/Elements/Penta.h =================================================================== --- /issm/trunk/src/c/objects/Elements/Penta.h (revision 8375) +++ /issm/trunk/src/c/objects/Elements/Penta.h (revision 8376) @@ -129,5 +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,Parameters* parameters){_error_("not supported yet!");} + void UpdateForcing(int index, IoModel* iomodel,int step,double time, int FieldEnum,Parameters* parameters); 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 8375) +++ /issm/trunk/src/c/objects/Elements/Tria.cpp (revision 8376) @@ -2728,5 +2728,4 @@ /*get inputs configured too: */ this->inputs->Configure(parameters); - } Index: /issm/trunk/src/c/objects/Inputs/PentaVertexForcing.cpp =================================================================== --- /issm/trunk/src/c/objects/Inputs/PentaVertexForcing.cpp (revision 8376) +++ /issm/trunk/src/c/objects/Inputs/PentaVertexForcing.cpp (revision 8376) @@ -0,0 +1,498 @@ +/*!\file PentaVertexForcing.c + * \brief: implementation of the PentaVertexForcing 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" + +/*PentaVertexForcing constructors and destructor*/ +/*FUNCTION PentaVertexForcing::PentaVertexForcing(){{{1*/ +PentaVertexForcing::PentaVertexForcing(){ + enum_type=NoneEnum; + values=NULL; + numtimesteps=0; + parameters=NULL; + return; +} +/*}}}*/ +/*FUNCTION PentaVertexForcing::PentaVertexForcing(int in_enum_type,double* time0values,double time0,int numtimesteps,Parameters* parameters){{{1*/ +PentaVertexForcing::PentaVertexForcing(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*6*sizeof(double)); + this->timesteps=(double*)xmalloc(this->numtimesteps*sizeof(double)); + for(int i=0;i<6;i++) this->values[i]=time0values[i]; + this->timesteps[0]=time0; + + this->parameters=parameters; + +} +/*}}}*/ +/*FUNCTION PentaVertexForcing::PentaVertexForcing(int in_enum_type,double* in_values,double* in_time,int numtimesteps,Parameters* parameters){{{1*/ +PentaVertexForcing::PentaVertexForcing(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*6*sizeof(double)); + this->timesteps=(double*)xmalloc(this->numtimesteps*sizeof(double)); + + memcpy(this->values,in_values,numtimesteps*6*sizeof(double)); + memcpy(this->timesteps,in_time,numtimesteps*sizeof(double)); + + this->parameters=parameters; + +} +/*}}}*/ +/*FUNCTION PentaVertexForcing::~PentaVertexForcing(){{{1*/ +PentaVertexForcing::~PentaVertexForcing(){ + xfree((void**)&this->values); + xfree((void**)&this->timesteps); + parameters=NULL; + return; +} +/*}}}*/ +/*FUNCTION PentaVertexForcing::AddTimeValues(double* values,int step,double time);{{{1*/ +void PentaVertexForcing::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; + for(int i=0;i<6;i++) this->values[6*step+i]=values[i]; +} +/*}}}*/ + +/*Object virtual functions definitions:*/ +/*FUNCTION PentaVertexForcing::Echo {{{1*/ +void PentaVertexForcing::Echo(void){ + this->DeepEcho(); +} +/*}}}*/ +/*FUNCTION PentaVertexForcing::DeepEcho{{{1*/ +void PentaVertexForcing::DeepEcho(void){ + + int i; + + printf("PentaVertexForcing:\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 %g %g %g]\n",this->timesteps[i],this->values[6*i+0],this->values[6*i+1],this->values[6*i+2],this->values[6*i+3],this->values[6*i+4],this->values[6*i+5]); + } +} +/*}}}*/ +/*FUNCTION PentaVertexForcing::Id{{{1*/ +int PentaVertexForcing::Id(void){ return -1; } +/*}}}*/ +/*FUNCTION PentaVertexForcing::MyRank{{{1*/ +int PentaVertexForcing::MyRank(void){ + extern int my_rank; + return my_rank; +} +/*}}}*/ +/*FUNCTION PentaVertexForcing::Marshall{{{1*/ +void PentaVertexForcing::Marshall(char** pmarshalled_dataset){ + + char* marshalled_dataset=NULL; + int enum_value=0; + + /*recover marshalled_dataset: */ + marshalled_dataset=*pmarshalled_dataset; + + /*get enum value of PentaVertexForcing: */ + enum_value=PentaVertexForcingEnum; + + /*marshall enum: */ + memcpy(marshalled_dataset,&enum_value,sizeof(enum_value));marshalled_dataset+=sizeof(enum_value); + + /*marshall PentaVertexForcing 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*6*sizeof(double));marshalled_dataset+=numtimesteps*6*sizeof(double); + memcpy(marshalled_dataset,timesteps,numtimesteps*sizeof(double));marshalled_dataset+=numtimesteps*sizeof(double); + + *pmarshalled_dataset=marshalled_dataset; +} +/*}}}*/ +/*FUNCTION PentaVertexForcing::MarshallSize{{{1*/ +int PentaVertexForcing::MarshallSize(){ + + return + +sizeof(enum_type)+ + +sizeof(numtimesteps)+ + 6*numtimesteps*sizeof(double)+ + numtimesteps*sizeof(double)+ + +sizeof(int); //sizeof(int) for enum value +} +/*}}}*/ +/*FUNCTION PentaVertexForcing::Demarshall{{{1*/ +void PentaVertexForcing::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*6*sizeof(double)); + + memcpy(values,marshalled_dataset,6*numtimesteps*sizeof(double));marshalled_dataset+=6*numtimesteps*sizeof(double); + memcpy(timesteps,marshalled_dataset,numtimesteps*sizeof(double));marshalled_dataset+=numtimesteps*sizeof(double); + + /*return: */ + *pmarshalled_dataset=marshalled_dataset; + return; +} +/*}}}*/ +/*FUNCTION PentaVertexForcing::Enum{{{1*/ +int PentaVertexForcing::Enum(void){ + + return PentaVertexForcingEnum; + +} +/*}}}*/ +/*FUNCTION PentaVertexForcing::copy{{{1*/ +Object* PentaVertexForcing::copy() { + + return new PentaVertexForcing(this->enum_type,this->values,this->timesteps,this->numtimesteps,this->parameters); + +} +/*}}}*/ + +/*PentaVertexForcing management*/ +/*FUNCTION PentaVertexForcing::EnumType{{{1*/ +int PentaVertexForcing::EnumType(void){ + + return this->enum_type; + +} +/*}}}*/ +/*FUNCTION PentaVertexForcing::SpawnTriaInput{{{1*/ +Input* PentaVertexForcing::SpawnTriaInput(int* indices){ + + /*output*/ + TriaVertexForcing* outforcing=NULL; + double* triavalues=NULL; + int i,j; + + /*Build tria values*/ + triavalues=(double*)xmalloc(this->numtimesteps*3*sizeof(double)); + for(i=0;i<3;i++)for(j=0;jnumtimesteps;j++) triavalues[3*j+i]=this->values[6*j+indices[i]]; + + /*Create new Tria forcing (copy of current forcing)*/ + outforcing=new TriaVertexForcing(this->enum_type,triavalues,this->timesteps,this->numtimesteps,this->parameters); + + /*Assign output*/ + return outforcing; +} +/*}}}*/ +/*FUNCTION PentaVertexForcing::SpawnResult{{{1*/ +ElementResult* PentaVertexForcing::SpawnResult(int step, double time){ + + double pentavalues[6]; + + /*Ok, we want to spawn a PentaVertexElementResult. We have the time, just get + *the correct values at the three nodes: */ + this->GetTimeValues(&pentavalues[0],time); + + return new PentaVertexElementResult(this->enum_type,&pentavalues[0],step,time); + +} +/*}}}*/ + +/*Object functions*/ +/*FUNCTION PentaVertexForcing::GetParameterValue(double* pvalue,GaussTria* gauss){{{1*/ +void PentaVertexForcing::GetParameterValue(double* pvalue,GaussTria* gauss){ + + double time; + double pentavalues[6]; + + /*First, recover current time from parameters: */ + parameters->FindParam(&time,TimeEnum); + + /*Retrieve interpolated values for this time step: */ + this->GetTimeValues(&pentavalues[0],time); + + /*Call TriaRef function*/ + TriaRef::GetParameterValue(pvalue,&pentavalues[0],gauss); +} +/*}}}*/ +/*FUNCTION PentaVertexForcing::GetParameterDerivativeValue(double* p, double* xyz_list, GaussTria* gauss){{{1*/ +void PentaVertexForcing::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 PentaVertexForcing::ChangeEnum{{{1*/ +void PentaVertexForcing::ChangeEnum(int newenumtype){ + this->enum_type=newenumtype; +} +/*}}}*/ +/*FUNCTION PentaVertexForcing::GetParameterAverage{{{1*/ +void PentaVertexForcing::GetParameterAverage(double* pvalue){ + + double time; + double pentavalues[6]; + + /*First, recover current time from parameters: */ + parameters->FindParam(&time,TimeEnum); + + /*Retrieve interpolated values for this time step: */ + this->GetTimeValues(&pentavalues[0],time); + + *pvalue=1./6.*(pentavalues[0]+pentavalues[1]+pentavalues[2]+pentavalues[3]+pentavalues[4]+pentavalues[5]); +} +/*}}}*/ + +/*Intermediary*/ +/*FUNCTION PentaVertexForcing::SquareMin{{{1*/ +void PentaVertexForcing::SquareMin(double* psquaremin, bool process_units,Parameters* parameters){ + + int i; + const int numnodes=6; + 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(&pentavalues[0],time); + + for(int i=0;inorm) norm=fabs(pentavalues[i]); + return norm; +} +/*}}}*/ +/*FUNCTION PentaVertexForcing::Max{{{1*/ +double PentaVertexForcing::Max(void){ + + const int numnodes=6; + double max; + double time; + double pentavalues[6]; + + /*First, recover current time from parameters: */ + parameters->FindParam(&time,TimeEnum); + + /*Retrieve interpolated values for this time step: */ + this->GetTimeValues(&pentavalues[0],time); + + max=pentavalues[0]; + + for(int i=1;imax) max=pentavalues[i]; + } + return max; +} +/*}}}*/ +/*FUNCTION PentaVertexForcing::MaxAbs{{{1*/ +double PentaVertexForcing::MaxAbs(void){ + + const int numnodes=6; + double max; + double time; + double pentavalues[6]; + + /*First, recover current time from parameters: */ + parameters->FindParam(&time,TimeEnum); + + /*Retrieve interpolated values for this time step: */ + this->GetTimeValues(&pentavalues[0],time); + + + max=fabs(pentavalues[0]); + + for(int i=1;imax) max=fabs(pentavalues[i]); + } + return max; +} +/*}}}*/ +/*FUNCTION PentaVertexForcing::Min{{{1*/ +double PentaVertexForcing::Min(void){ + + const int numnodes=6; + double min; + double time; + double pentavalues[6]; + + /*First, recover current time from parameters: */ + parameters->FindParam(&time,TimeEnum); + + /*Retrieve interpolated values for this time step: */ + this->GetTimeValues(&pentavalues[0],time); + + min=pentavalues[0]; + + for(int i=1;iFindParam(&time,TimeEnum); + + /*Retrieve interpolated values for this time step: */ + this->GetTimeValues(&pentavalues[0],time); + + min=fabs(pentavalues[0]); + + for(int i=1;iFindParam(&time,TimeEnum); + + /*Retrieve interpolated values for this time step: */ + this->GetTimeValues(&pentavalues[0],time); + + VecSetValues(vector,numvertices,doflist,(const double*)&pentavalues[0],INSERT_VALUES); + +} /*}}}*/ +/*FUNCTION PentaVertexForcing::GetTimeValues{{{1*/ +void PentaVertexForcing::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: */ + for(i=0;i<6;i++) outvalues[i]=this->values[6*0+i]; + found=true; + } + else if(intime>this->timesteps[this->numtimesteps-1]){ + /*get values for the last time: */ + for(i=0;i<6;i++) outvalues[i]=this->values[6*(this->numtimesteps-1)+i]; + 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: */ + for(j=0;j<6;j++) outvalues[j]=this->values[6*i+j]; + 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<6;j++) outvalues[j]=this->values[6*i+j]+alpha*(this->values[6*(i+1)+j]-this->values[6*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 PentaVertexForcing::Configure{{{1*/ +void PentaVertexForcing::Configure(Parameters* parameters){ + printf("Configuring %p\n",parameters); + this->parameters=parameters; +} +/*}}}*/ Index: /issm/trunk/src/c/objects/Inputs/PentaVertexForcing.h =================================================================== --- /issm/trunk/src/c/objects/Inputs/PentaVertexForcing.h (revision 8376) +++ /issm/trunk/src/c/objects/Inputs/PentaVertexForcing.h (revision 8376) @@ -0,0 +1,91 @@ +/*! \file PentaVertexForcing.h + * \brief: header file for triavertexforcing object + */ + + +#ifndef _PENTAVERTEXFORCING_H_ +#define _PENTAVERTEXFORCING_H_ + +/*Headers:*/ +/*{{{1*/ +#include "./Input.h" +#include "../Elements/TriaRef.h" +class GaussTria; +class Parameters; +/*}}}*/ + +class PentaVertexForcing: public Input,public TriaRef{ + + public: + /*just hold 6 values for 6 vertices, for all time steps: */ + int enum_type; + int numtimesteps; + double* values; + double* timesteps; + Parameters* parameters; //to find current time. + + /*PentaVertexForcing constructors, destructors: {{{1*/ + PentaVertexForcing(); + PentaVertexForcing(int enum_type,double* valuestime0,double time0,int numtimesteps,Parameters* parameters); + PentaVertexForcing(int enum_type,double* values,double* timesteps,int numtimesteps,Parameters* parameters); + ~PentaVertexForcing(); + 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(); + /*}}}*/ + /*PentaVertexForcing 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/objects.h =================================================================== --- /issm/trunk/src/c/objects/objects.h (revision 8375) +++ /issm/trunk/src/c/objects/objects.h (revision 8376) @@ -79,4 +79,5 @@ #include "./Inputs/IntInput.h" #include "./Inputs/PentaVertexInput.h" +#include "./Inputs/PentaVertexForcing.h" #include "./Inputs/TriaVertexInput.h" #include "./Inputs/TriaVertexForcing.h"