Index: /issm/trunk/src/ad/Makefile.am =================================================================== --- /issm/trunk/src/ad/Makefile.am (revision 9816) +++ /issm/trunk/src/ad/Makefile.am (revision 9817) @@ -260,6 +260,4 @@ ./modules/NodesDofx/NodesDofx.h\ ./modules/NodesDofx/NodesDofx.cpp\ - ./modules/NodalValuex/NodalValuex.h\ - ./modules/NodalValuex/NodalValuex.cpp\ ./modules/VerticesDofx/VerticesDofx.h\ ./modules/VerticesDofx/VerticesDofx.cpp\ @@ -533,4 +531,6 @@ ./modules/MassFluxx/MassFluxx.cpp\ ./modules/MassFluxx/MassFluxx.h\ + ./modules/NodalValuex/NodalValuex.h\ + ./modules/NodalValuex/NodalValuex.cpp\ ./modules/SurfaceAreax/SurfaceAreax.h\ ./modules/SurfaceAreax/SurfaceAreax.cpp\ Index: /issm/trunk/src/c/objects/Elements/Element.h =================================================================== --- /issm/trunk/src/c/objects/Elements/Element.h (revision 9816) +++ /issm/trunk/src/c/objects/Elements/Element.h (revision 9817) @@ -40,4 +40,53 @@ virtual void GetParameterListOnVertices(double* pvalue,int enumtype,double defaultvalue)=0; virtual void GetParameterValue(double* pvalue,Node* node,int enumtype)=0; + + virtual double SurfaceArea(void)=0; + virtual void InputDepthAverageAtBase(int enum_type,int average_enum_type,int object_enum)=0; + virtual void ComputeBasalStress(Vec sigma_b)=0; + virtual void ComputeStrainRate(Vec eps)=0; + virtual void PatchSize(int* pnumrows, int* pnumvertices,int* pnumnodes)=0; + virtual void PatchFill(int* pcount, Patch* patch)=0; + virtual void DeleteResults(void)=0; + virtual void Update(int index, IoModel* iomodel,int analysis_counter,int analysis_type)=0; + virtual void InputToResult(int enum_type,int step,double time)=0; + virtual void InputDuplicate(int original_enum,int new_enum)=0; + virtual void InputCreate(double scalar,int name,int code)=0; + virtual void InputCreate(double* vector, int index,IoModel* iomodel,int M,int N,int vector_type,int vector_enum,int code)=0; + virtual void ProcessResultsUnits(void)=0; + virtual void RequestedOutput(int output_enum,int step,double time)=0; + + virtual void InputScale(int enum_type,double scale_factor)=0; + virtual void GetVectorFromInputs(Vec vector,int NameEnum)=0; + virtual void InputArtificialNoise(int enum_type,double min,double max)=0; + virtual bool InputConvergence(double* eps, int* enums,int num_enums,int* criterionenums,double* criterionvalues,int num_criterionenums)=0; + virtual void AverageOntoPartition(Vec partition_contributions,Vec partition_areas,double* vertex_response,double* qmu_part)=0; + virtual int* GetHorizontalNeighboorSids(void)=0; + virtual double TimeAdapt()=0; + virtual void AgressiveMigration()=0; + virtual void SoftMigration(double* sheet_ungrounding)=0; + virtual void ShelfSync()=0; + virtual void PotentialSheetUngrounding(Vec potential_sheet_ungrounding)=0; + virtual void MigrateGroundingLine()=0; + virtual int UpdateShelfStatus(Vec new_shelf_nodes)=0; + virtual void UpdateShelfFlags(double* new_shelf_nodes)=0; + virtual int UpdatePotentialSheetUngrounding(double* potential_sheet_ungrounding,Vec vec_nodes_on_iceshelf,double* nodes_on_iceshelf)=0; + + #ifdef _HAVE_RESPONSES_ + virtual void MinVel(double* pminvel, bool process_units)=0; + virtual void MaxVel(double* pmaxvel, bool process_units)=0; + virtual void MinVx(double* pminvx, bool process_units)=0; + virtual void MaxVx(double* pmaxvx, bool process_units)=0; + virtual void MaxAbsVx(double* pmaxabsvx, bool process_units)=0; + virtual void MinVy(double* pminvy, bool process_units)=0; + virtual void MaxVy(double* pmaxvy, bool process_units)=0; + virtual void MaxAbsVy(double* pmaxabsvy, bool process_units)=0; + virtual void MinVz(double* pminvz, bool process_units)=0; + virtual void MaxVz(double* pmaxvz, bool process_units)=0; + virtual void MaxAbsVz(double* pmaxabsvz, bool process_units)=0; + virtual double MassFlux(double* segment,bool process_units)=0; + virtual double RheologyBbarx(void)=0; + virtual int NodalValue(double* pvalue, int index, int natureofdataenum,bool process_units)=0; + #endif + #ifdef _HAVE_CONTROL_ virtual void Gradj(Vec gradient,int control_type)=0; @@ -56,50 +105,4 @@ virtual void InputControlUpdate(double scalar,bool save_parameter)=0; #endif - virtual double SurfaceArea(void)=0; - virtual void InputDepthAverageAtBase(int enum_type,int average_enum_type,int object_enum)=0; - virtual void ComputeBasalStress(Vec sigma_b)=0; - virtual void ComputeStrainRate(Vec eps)=0; - virtual double MassFlux(double* segment,bool process_units)=0; - virtual void PatchSize(int* pnumrows, int* pnumvertices,int* pnumnodes)=0; - virtual void PatchFill(int* pcount, Patch* patch)=0; - virtual void DeleteResults(void)=0; - virtual void Update(int index, IoModel* iomodel,int analysis_counter,int analysis_type)=0; - virtual void InputToResult(int enum_type,int step,double time)=0; - virtual void InputDuplicate(int original_enum,int new_enum)=0; - virtual void InputCreate(double scalar,int name,int code)=0; - virtual void InputCreate(double* vector, int index,IoModel* iomodel,int M,int N,int vector_type,int vector_enum,int code)=0; - virtual void ProcessResultsUnits(void)=0; - virtual void RequestedOutput(int output_enum,int step,double time)=0; - virtual void MinVel(double* pminvel, bool process_units)=0; - virtual void MaxVel(double* pmaxvel, bool process_units)=0; - virtual void MinVx(double* pminvx, bool process_units)=0; - virtual void MaxVx(double* pmaxvx, bool process_units)=0; - virtual void MaxAbsVx(double* pmaxabsvx, bool process_units)=0; - virtual void MinVy(double* pminvy, bool process_units)=0; - virtual void MaxVy(double* pmaxvy, bool process_units)=0; - virtual void MaxAbsVy(double* pmaxabsvy, bool process_units)=0; - virtual void MinVz(double* pminvz, bool process_units)=0; - virtual void MaxVz(double* pmaxvz, bool process_units)=0; - virtual void MaxAbsVz(double* pmaxabsvz, bool process_units)=0; - virtual void InputScale(int enum_type,double scale_factor)=0; - virtual void GetVectorFromInputs(Vec vector,int NameEnum)=0; - virtual void InputArtificialNoise(int enum_type,double min,double max)=0; - virtual bool InputConvergence(double* eps, int* enums,int num_enums,int* criterionenums,double* criterionvalues,int num_criterionenums)=0; - virtual void AverageOntoPartition(Vec partition_contributions,Vec partition_areas,double* vertex_response,double* qmu_part)=0; - virtual int* GetHorizontalNeighboorSids(void)=0; - virtual double TimeAdapt()=0; - virtual void AgressiveMigration()=0; - virtual void SoftMigration(double* sheet_ungrounding)=0; - virtual void ShelfSync()=0; - virtual void PotentialSheetUngrounding(Vec potential_sheet_ungrounding)=0; - virtual void MigrateGroundingLine()=0; - virtual int UpdateShelfStatus(Vec new_shelf_nodes)=0; - virtual void UpdateShelfFlags(double* new_shelf_nodes)=0; - virtual int UpdatePotentialSheetUngrounding(double* potential_sheet_ungrounding,Vec vec_nodes_on_iceshelf,double* nodes_on_iceshelf)=0; - virtual double RheologyBbarx(void)=0; - virtual int NodalValue(double* pvalue, int index, int natureofdataenum,bool process_units)=0; - - /*Implementation: */ - }; #endif Index: /issm/trunk/src/c/objects/Elements/Penta.cpp =================================================================== --- /issm/trunk/src/c/objects/Elements/Penta.cpp (revision 9816) +++ /issm/trunk/src/c/objects/Elements/Penta.cpp (revision 9817) @@ -2001,101 +2001,4 @@ } /*}}}*/ -/*FUNCTION Penta::MassFlux {{{1*/ -double Penta::MassFlux( double* segment,bool process_units){ - _error_(" not supported yet!"); -} -/*}}}*/ -/*FUNCTION Penta::MaxAbsVx{{{1*/ -void Penta::MaxAbsVx(double* pmaxabsvx, bool process_units){ - - /*Get maximum:*/ - double maxabsvx=this->inputs->MaxAbs(VxEnum); - - /*process units if requested: */ - if(process_units) maxabsvx=UnitConversion(maxabsvx,IuToExtEnum,VxEnum); - - /*Assign output pointers:*/ - *pmaxabsvx=maxabsvx; -} -/*}}}*/ -/*FUNCTION Penta::MaxAbsVy{{{1*/ -void Penta::MaxAbsVy(double* pmaxabsvy, bool process_units){ - - /*Get maximum:*/ - double maxabsvy=this->inputs->MaxAbs(VyEnum); - - /*process units if requested: */ - if(process_units) maxabsvy=UnitConversion(maxabsvy,IuToExtEnum,VyEnum); - - /*Assign output pointers:*/ - *pmaxabsvy=maxabsvy; -} -/*}}}*/ -/*FUNCTION Penta::MaxAbsVz{{{1*/ -void Penta::MaxAbsVz(double* pmaxabsvz, bool process_units){ - - /*Get maximum:*/ - double maxabsvz=this->inputs->MaxAbs(VzEnum); - - /*process units if requested: */ - if(process_units) maxabsvz=UnitConversion(maxabsvz,IuToExtEnum,VyEnum); - - /*Assign output pointers:*/ - *pmaxabsvz=maxabsvz; -} -/*}}}*/ -/*FUNCTION Penta::MaxVel{{{1*/ -void Penta::MaxVel(double* pmaxvel, bool process_units){ - - /*Get maximum:*/ - double maxvel=this->inputs->Max(VelEnum); - - /*process units if requested: */ - if(process_units) maxvel=UnitConversion(maxvel,IuToExtEnum,VelEnum); - - /*Assign output pointers:*/ - *pmaxvel=maxvel; - -} -/*}}}*/ -/*FUNCTION Penta::MaxVx{{{1*/ -void Penta::MaxVx(double* pmaxvx, bool process_units){ - - /*Get maximum:*/ - double maxvx=this->inputs->Max(VxEnum); - - /*process units if requested: */ - if(process_units) maxvx=UnitConversion(maxvx,IuToExtEnum,VxEnum); - - /*Assign output pointers:*/ - *pmaxvx=maxvx; -} -/*}}}*/ -/*FUNCTION Penta::MaxVy{{{1*/ -void Penta::MaxVy(double* pmaxvy, bool process_units){ - - /*Get maximum:*/ - double maxvy=this->inputs->Max(VyEnum); - - /*process units if requested: */ - if(process_units) maxvy=UnitConversion(maxvy,IuToExtEnum,VyEnum); - - /*Assign output pointers:*/ - *pmaxvy=maxvy; -} -/*}}}*/ -/*FUNCTION Penta::MaxVz{{{1*/ -void Penta::MaxVz(double* pmaxvz, bool process_units){ - - /*Get maximum:*/ - double maxvz=this->inputs->Max(VzEnum); - - /*process units if requested: */ - if(process_units) maxvz=UnitConversion(maxvz,IuToExtEnum,VzEnum); - - /*Assign output pointers:*/ - *pmaxvz=maxvz; -} -/*}}}*/ /*FUNCTION Penta::MigrateGroundingLine{{{1*/ void Penta::MigrateGroundingLine(void){ @@ -2125,91 +2028,8 @@ } /*}}}*/ -/*FUNCTION Penta::MinVel{{{1*/ -void Penta::MinVel(double* pminvel, bool process_units){ - - /*Get minimum:*/ - double minvel=this->inputs->Min(VelEnum); - - /*process units if requested: */ - if(process_units) minvel=UnitConversion(minvel,IuToExtEnum,VelEnum); - - /*Assign output pointers:*/ - *pminvel=minvel; -} -/*}}}*/ -/*FUNCTION Penta::MinVx{{{1*/ -void Penta::MinVx(double* pminvx, bool process_units){ - - /*Get minimum:*/ - double minvx=this->inputs->Min(VxEnum); - - /*process units if requested: */ - if(process_units) minvx=UnitConversion(minvx,IuToExtEnum,VxEnum); - - /*Assign output pointers:*/ - *pminvx=minvx; -} -/*}}}*/ -/*FUNCTION Penta::MinVy{{{1*/ -void Penta::MinVy(double* pminvy, bool process_units){ - - /*Get minimum:*/ - double minvy=this->inputs->Min(VyEnum); - - /*process units if requested: */ - if(process_units) minvy=UnitConversion(minvy,IuToExtEnum,VyEnum); - - /*Assign output pointers:*/ - *pminvy=minvy; -} -/*}}}*/ -/*FUNCTION Penta::MinVz{{{1*/ -void Penta::MinVz(double* pminvz, bool process_units){ - - /*Get minimum:*/ - double minvz=this->inputs->Min(VzEnum); - - /*process units if requested: */ - if(process_units) minvz=UnitConversion(minvz,IuToExtEnum,VzEnum); - - /*Assign output pointers:*/ - *pminvz=minvz; -} -/*}}}*/ /*FUNCTION Penta::MyRank {{{1*/ int Penta::MyRank(void){ extern int my_rank; return my_rank; -} -/*}}}*/ -/*FUNCTION Penta::NodalValue {{{1*/ -int Penta::NodalValue(double* pvalue, int index, int natureofdataenum,bool process_units){ - - int i; - int found=0; - double value; - Input* data=NULL; - GaussPenta* gauss=NULL; - - /*First, serarch the input: */ - data=inputs->GetInput(natureofdataenum); - - /*figure out if we have the vertex id: */ - found=0; - for(i=0;iGetVertexId()){ - /*Do we have natureofdataenum in our inputs? :*/ - if(data){ - /*ok, we are good. retrieve value of input at vertex :*/ - gauss=new GaussPenta(); gauss->GaussVertex(i); - data->GetParameterValue(&value,gauss); - found=1; - break; - } - } - } - - if(found)*pvalue=value; - return found; } /*}}}*/ @@ -2397,32 +2217,4 @@ } -} -/*}}}*/ -/*FUNCTION Penta::RheologyBbarx{{{1*/ -double Penta::RheologyBbarx(void){ - - double Bbar; - Penta* base=NULL; - - /*We here compute the vertically averaged Bbar. To do so, we go to the - * base element under this one, and call this->InputDepthAverageAtBase:*/ - - base=this; - - if(!base->IsOnBed()){ - - for(;;){ - base=base->verticalneighbors[0]; - if(base->IsOnBed())break; - } - } - - /*We have found the base under this. Call InputDepthAverageAtBase on it: */ - base->InputDepthAverageAtBase(MaterialsRheologyBEnum,MaterialsRheologyBbarEnum,MaterialsEnum); - - /*Ok, base->matice now has the Bbar. Ask Bbar from this enum: */ - Bbar=base->matice->GetBbar(); - - return Bbar; } /*}}}*/ @@ -2821,4 +2613,215 @@ /*}}}*/ +#ifdef _HAVE_RESPONSES_ +/*FUNCTION Penta::RheologyBbarx{{{1*/ +double Penta::RheologyBbarx(void){ + + double Bbar; + Penta* base=NULL; + + /*We here compute the vertically averaged Bbar. To do so, we go to the + * base element under this one, and call this->InputDepthAverageAtBase:*/ + + base=this; + + if(!base->IsOnBed()){ + + for(;;){ + base=base->verticalneighbors[0]; + if(base->IsOnBed())break; + } + } + + /*We have found the base under this. Call InputDepthAverageAtBase on it: */ + base->InputDepthAverageAtBase(MaterialsRheologyBEnum,MaterialsRheologyBbarEnum,MaterialsEnum); + + /*Ok, base->matice now has the Bbar. Ask Bbar from this enum: */ + Bbar=base->matice->GetBbar(); + + return Bbar; +} +/*}}}*/ +/*FUNCTION Penta::NodalValue {{{1*/ +int Penta::NodalValue(double* pvalue, int index, int natureofdataenum,bool process_units){ + + int i; + int found=0; + double value; + Input* data=NULL; + GaussPenta* gauss=NULL; + + /*First, serarch the input: */ + data=inputs->GetInput(natureofdataenum); + + /*figure out if we have the vertex id: */ + found=0; + for(i=0;iGetVertexId()){ + /*Do we have natureofdataenum in our inputs? :*/ + if(data){ + /*ok, we are good. retrieve value of input at vertex :*/ + gauss=new GaussPenta(); gauss->GaussVertex(i); + data->GetParameterValue(&value,gauss); + found=1; + break; + } + } + } + + if(found)*pvalue=value; + return found; +} +/*}}}*/ +/*FUNCTION Penta::MinVel{{{1*/ +void Penta::MinVel(double* pminvel, bool process_units){ + + /*Get minimum:*/ + double minvel=this->inputs->Min(VelEnum); + + /*process units if requested: */ + if(process_units) minvel=UnitConversion(minvel,IuToExtEnum,VelEnum); + + /*Assign output pointers:*/ + *pminvel=minvel; +} +/*}}}*/ +/*FUNCTION Penta::MinVx{{{1*/ +void Penta::MinVx(double* pminvx, bool process_units){ + + /*Get minimum:*/ + double minvx=this->inputs->Min(VxEnum); + + /*process units if requested: */ + if(process_units) minvx=UnitConversion(minvx,IuToExtEnum,VxEnum); + + /*Assign output pointers:*/ + *pminvx=minvx; +} +/*}}}*/ +/*FUNCTION Penta::MinVy{{{1*/ +void Penta::MinVy(double* pminvy, bool process_units){ + + /*Get minimum:*/ + double minvy=this->inputs->Min(VyEnum); + + /*process units if requested: */ + if(process_units) minvy=UnitConversion(minvy,IuToExtEnum,VyEnum); + + /*Assign output pointers:*/ + *pminvy=minvy; +} +/*}}}*/ +/*FUNCTION Penta::MinVz{{{1*/ +void Penta::MinVz(double* pminvz, bool process_units){ + + /*Get minimum:*/ + double minvz=this->inputs->Min(VzEnum); + + /*process units if requested: */ + if(process_units) minvz=UnitConversion(minvz,IuToExtEnum,VzEnum); + + /*Assign output pointers:*/ + *pminvz=minvz; +} +/*}}}*/ +/*FUNCTION Penta::MassFlux {{{1*/ +double Penta::MassFlux( double* segment,bool process_units){ + _error_(" not supported yet!"); +} +/*}}}*/ +/*FUNCTION Penta::MaxAbsVx{{{1*/ +void Penta::MaxAbsVx(double* pmaxabsvx, bool process_units){ + + /*Get maximum:*/ + double maxabsvx=this->inputs->MaxAbs(VxEnum); + + /*process units if requested: */ + if(process_units) maxabsvx=UnitConversion(maxabsvx,IuToExtEnum,VxEnum); + + /*Assign output pointers:*/ + *pmaxabsvx=maxabsvx; +} +/*}}}*/ +/*FUNCTION Penta::MaxAbsVy{{{1*/ +void Penta::MaxAbsVy(double* pmaxabsvy, bool process_units){ + + /*Get maximum:*/ + double maxabsvy=this->inputs->MaxAbs(VyEnum); + + /*process units if requested: */ + if(process_units) maxabsvy=UnitConversion(maxabsvy,IuToExtEnum,VyEnum); + + /*Assign output pointers:*/ + *pmaxabsvy=maxabsvy; +} +/*}}}*/ +/*FUNCTION Penta::MaxAbsVz{{{1*/ +void Penta::MaxAbsVz(double* pmaxabsvz, bool process_units){ + + /*Get maximum:*/ + double maxabsvz=this->inputs->MaxAbs(VzEnum); + + /*process units if requested: */ + if(process_units) maxabsvz=UnitConversion(maxabsvz,IuToExtEnum,VyEnum); + + /*Assign output pointers:*/ + *pmaxabsvz=maxabsvz; +} +/*}}}*/ +/*FUNCTION Penta::MaxVel{{{1*/ +void Penta::MaxVel(double* pmaxvel, bool process_units){ + + /*Get maximum:*/ + double maxvel=this->inputs->Max(VelEnum); + + /*process units if requested: */ + if(process_units) maxvel=UnitConversion(maxvel,IuToExtEnum,VelEnum); + + /*Assign output pointers:*/ + *pmaxvel=maxvel; + +} +/*}}}*/ +/*FUNCTION Penta::MaxVx{{{1*/ +void Penta::MaxVx(double* pmaxvx, bool process_units){ + + /*Get maximum:*/ + double maxvx=this->inputs->Max(VxEnum); + + /*process units if requested: */ + if(process_units) maxvx=UnitConversion(maxvx,IuToExtEnum,VxEnum); + + /*Assign output pointers:*/ + *pmaxvx=maxvx; +} +/*}}}*/ +/*FUNCTION Penta::MaxVy{{{1*/ +void Penta::MaxVy(double* pmaxvy, bool process_units){ + + /*Get maximum:*/ + double maxvy=this->inputs->Max(VyEnum); + + /*process units if requested: */ + if(process_units) maxvy=UnitConversion(maxvy,IuToExtEnum,VyEnum); + + /*Assign output pointers:*/ + *pmaxvy=maxvy; +} +/*}}}*/ +/*FUNCTION Penta::MaxVz{{{1*/ +void Penta::MaxVz(double* pmaxvz, bool process_units){ + + /*Get maximum:*/ + double maxvz=this->inputs->Max(VzEnum); + + /*process units if requested: */ + if(process_units) maxvz=UnitConversion(maxvz,IuToExtEnum,VzEnum); + + /*Assign output pointers:*/ + *pmaxvz=maxvz; +} +/*}}}*/ +#endif + #ifdef _HAVE_THERMAL_ /*FUNCTION Penta::CreateKMatrixEnthalpy {{{1*/ Index: /issm/trunk/src/c/objects/Elements/Penta.h =================================================================== --- /issm/trunk/src/c/objects/Elements/Penta.h (revision 9816) +++ /issm/trunk/src/c/objects/Elements/Penta.h (revision 9817) @@ -102,12 +102,4 @@ void InputToResult(int enum_type,int step,double time); - double MassFlux(double* segment,bool process_units); - void MaxAbsVx(double* pmaxabsvx, bool process_units); - void MaxAbsVy(double* pmaxabsvy, bool process_units); - void MaxAbsVz(double* pmaxabsvz, bool process_units); - void MaxVel(double* pmaxvel, bool process_units); - void MaxVx(double* pmaxvx, bool process_units); - void MaxVy(double* pmaxvy, bool process_units); - void MaxVz(double* pmaxvz, bool process_units); void AgressiveMigration(); void SoftMigration(double* sheet_ungrounding); @@ -116,8 +108,4 @@ void RequestedOutput(int output_enum,int step,double time); void MigrateGroundingLine(); - void MinVel(double* pminvel, bool process_units); - void MinVx(double* pminvx, bool process_units); - void MinVy(double* pminvy, bool process_units); - void MinVz(double* pminvz, bool process_units); void PatchFill(int* pcount, Patch* patch); void PatchSize(int* pnumrows, int* pnumvertices,int* pnumnodes); @@ -130,7 +118,22 @@ double TimeAdapt(); int* GetHorizontalNeighboorSids(void); + void ViscousHeatingCreateInput(void); + + #ifdef _HAVE_RESPONSES_ + void MinVel(double* pminvel, bool process_units); + void MinVx(double* pminvx, bool process_units); + void MinVy(double* pminvy, bool process_units); + void MinVz(double* pminvz, bool process_units); double RheologyBbarx(void); - void ViscousHeatingCreateInput(void); int NodalValue(double* pvalue, int index, int natureofdataenum,bool process_units); + double MassFlux(double* segment,bool process_units); + void MaxAbsVx(double* pmaxabsvx, bool process_units); + void MaxAbsVy(double* pmaxabsvy, bool process_units); + void MaxAbsVz(double* pmaxabsvz, bool process_units); + void MaxVel(double* pmaxvel, bool process_units); + void MaxVx(double* pmaxvx, bool process_units); + void MaxVy(double* pmaxvy, bool process_units); + void MaxVz(double* pmaxvz, bool process_units); + #endif #ifdef _HAVE_CONTROL_ Index: /issm/trunk/src/c/objects/Elements/Tria.cpp =================================================================== --- /issm/trunk/src/c/objects/Elements/Tria.cpp (revision 9816) +++ /issm/trunk/src/c/objects/Elements/Tria.cpp (revision 9817) @@ -1926,159 +1926,4 @@ } /*}}}*/ -/*FUNCTION Tria::MassFlux {{{1*/ -double Tria::MassFlux( double* segment,bool process_units){ - - const int numdofs=2; - - int i; - double mass_flux=0; - double xyz_list[NUMVERTICES][3]; - double normal[2]; - double length,rho_ice; - double x1,y1,x2,y2,h1,h2; - double vx1,vx2,vy1,vy2; - GaussTria* gauss_1=NULL; - GaussTria* gauss_2=NULL; - - /*Get material parameters :*/ - rho_ice=matpar->GetRhoIce(); - - /*First off, check that this segment belongs to this element: */ - if ((int)*(segment+4)!=this->id)_error_("%s%i%s%i","error message: segment with id ",(int)*(segment+4)," does not belong to element with id:",this->id); - - /*Recover segment node locations: */ - x1=*(segment+0); y1=*(segment+1); x2=*(segment+2); y2=*(segment+3); - - /*Get xyz list: */ - GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES); - - /*get area coordinates of 0 and 1 locations: */ - gauss_1=new GaussTria(); - gauss_1->GaussFromCoords(x1,y1,&xyz_list[0][0]); - gauss_2=new GaussTria(); - gauss_2->GaussFromCoords(x2,y2,&xyz_list[0][0]); - - normal[0]=cos(atan2(x1-x2,y2-y1)); - normal[1]=sin(atan2(x1-x2,y2-y1)); - - length=sqrt(pow(x2-x1,2.0)+pow(y2-y1,2)); - - Input* thickness_input=inputs->GetInput(ThicknessEnum); _assert_(thickness_input); - Input* vx_input=inputs->GetInput(VxEnum); _assert_(vx_input); - Input* vy_input=inputs->GetInput(VyEnum); _assert_(vy_input); - - thickness_input->GetParameterValue(&h1, gauss_1); - thickness_input->GetParameterValue(&h2, gauss_2); - vx_input->GetParameterValue(&vx1,gauss_1); - vx_input->GetParameterValue(&vx2,gauss_2); - vy_input->GetParameterValue(&vy1,gauss_1); - vy_input->GetParameterValue(&vy2,gauss_2); - - mass_flux= rho_ice*length*( - (ONETHIRD*(h1-h2)*(vx1-vx2)+0.5*h2*(vx1-vx2)+0.5*(h1-h2)*vx2+h2*vx2)*normal[0]+ - (ONETHIRD*(h1-h2)*(vy1-vy2)+0.5*h2*(vy1-vy2)+0.5*(h1-h2)*vy2+h2*vy2)*normal[1] - ); - - /*Process units: */ - mass_flux=UnitConversion(mass_flux,IuToExtEnum,MassFluxEnum); - - /*clean up and return:*/ - delete gauss_1; - delete gauss_2; - return mass_flux; -} -/*}}}*/ -/*FUNCTION Tria::MaxAbsVx{{{1*/ -void Tria::MaxAbsVx(double* pmaxabsvx, bool process_units){ - - /*Get maximum:*/ - double maxabsvx=this->inputs->MaxAbs(VxEnum); - - /*process units if requested: */ - if(process_units) maxabsvx=UnitConversion(maxabsvx,IuToExtEnum,VxEnum); - - /*Assign output pointers:*/ - *pmaxabsvx=maxabsvx; -} -/*}}}*/ -/*FUNCTION Tria::MaxAbsVy{{{1*/ -void Tria::MaxAbsVy(double* pmaxabsvy, bool process_units){ - - /*Get maximum:*/ - double maxabsvy=this->inputs->MaxAbs(VyEnum); - - /*process units if requested: */ - if(process_units) maxabsvy=UnitConversion(maxabsvy,IuToExtEnum,VyEnum); - - /*Assign output pointers:*/ - *pmaxabsvy=maxabsvy; -} -/*}}}*/ -/*FUNCTION Tria::MaxAbsVz{{{1*/ -void Tria::MaxAbsVz(double* pmaxabsvz, bool process_units){ - - /*Get maximum:*/ - double maxabsvz=this->inputs->MaxAbs(VzEnum); - - /*process units if requested: */ - if(process_units) maxabsvz=UnitConversion(maxabsvz,IuToExtEnum,VyEnum); - - /*Assign output pointers:*/ - *pmaxabsvz=maxabsvz; -} -/*}}}*/ -/*FUNCTION Tria::MaxVel{{{1*/ -void Tria::MaxVel(double* pmaxvel, bool process_units){ - - /*Get maximum:*/ - double maxvel=this->inputs->Max(VelEnum); - - /*process units if requested: */ - if(process_units) maxvel=UnitConversion(maxvel,IuToExtEnum,VelEnum); - - /*Assign output pointers:*/ - *pmaxvel=maxvel; -} -/*}}}*/ -/*FUNCTION Tria::MaxVx{{{1*/ -void Tria::MaxVx(double* pmaxvx, bool process_units){ - - /*Get maximum:*/ - double maxvx=this->inputs->Max(VxEnum); - - /*process units if requested: */ - if(process_units) maxvx=UnitConversion(maxvx,IuToExtEnum,VxEnum); - - /*Assign output pointers:*/ - *pmaxvx=maxvx; -} -/*}}}*/ -/*FUNCTION Tria::MaxVy{{{1*/ -void Tria::MaxVy(double* pmaxvy, bool process_units){ - - /*Get maximum:*/ - double maxvy=this->inputs->Max(VyEnum); - - /*process units if requested: */ - if(process_units) maxvy=UnitConversion(maxvy,IuToExtEnum,VyEnum); - - /*Assign output pointers:*/ - *pmaxvy=maxvy; - -} -/*}}}*/ -/*FUNCTION Tria::MaxVz{{{1*/ -void Tria::MaxVz(double* pmaxvz, bool process_units){ - - /*Get maximum:*/ - double maxvz=this->inputs->Max(VzEnum); - - /*process units if requested: */ - if(process_units) maxvz=UnitConversion(maxvz,IuToExtEnum,VzEnum); - - /*Assign output pointers:*/ - *pmaxvz=maxvz; -} -/*}}}*/ /*FUNCTION Tria::MigrateGroundingLine{{{1*/ void Tria::MigrateGroundingLine(void){ @@ -2156,91 +2001,8 @@ } /*}}}*/ -/*FUNCTION Tria::MinVel{{{1*/ -void Tria::MinVel(double* pminvel, bool process_units){ - - /*Get minimum:*/ - double minvel=this->inputs->Min(VelEnum); - - /*process units if requested: */ - if(process_units) minvel=UnitConversion(minvel,IuToExtEnum,VelEnum); - - /*Assign output pointers:*/ - *pminvel=minvel; -} -/*}}}*/ -/*FUNCTION Tria::MinVx{{{1*/ -void Tria::MinVx(double* pminvx, bool process_units){ - - /*Get minimum:*/ - double minvx=this->inputs->Min(VxEnum); - - /*process units if requested: */ - if(process_units) minvx=UnitConversion(minvx,IuToExtEnum,VxEnum); - - /*Assign output pointers:*/ - *pminvx=minvx; -} -/*}}}*/ -/*FUNCTION Tria::MinVy{{{1*/ -void Tria::MinVy(double* pminvy, bool process_units){ - - /*Get minimum:*/ - double minvy=this->inputs->Min(VyEnum); - - /*process units if requested: */ - if(process_units) minvy=UnitConversion(minvy,IuToExtEnum,VyEnum); - - /*Assign output pointers:*/ - *pminvy=minvy; -} -/*}}}*/ -/*FUNCTION Tria::MinVz{{{1*/ -void Tria::MinVz(double* pminvz, bool process_units){ - - /*Get minimum:*/ - double minvz=this->inputs->Min(VzEnum); - - /*process units if requested: */ - if(process_units) minvz=UnitConversion(minvz,IuToExtEnum,VzEnum); - - /*Assign output pointers:*/ - *pminvz=minvz; -} -/*}}}*/ /*FUNCTION Tria::MyRank {{{1*/ int Tria::MyRank(void){ extern int my_rank; return my_rank; -} -/*}}}*/ -/*FUNCTION Tria::NodalValue {{{1*/ -int Tria::NodalValue(double* pvalue, int index, int natureofdataenum,bool process_units){ - - int i; - int found=0; - double value; - Input* data=NULL; - GaussTria *gauss = NULL; - - /*First, serarch the input: */ - data=inputs->GetInput(natureofdataenum); - - /*figure out if we have the vertex id: */ - found=0; - for(i=0;iGetVertexId()){ - /*Do we have natureofdataenum in our inputs? :*/ - if(data){ - /*ok, we are good. retrieve value of input at vertex :*/ - gauss=new GaussTria(); gauss->GaussVertex(i); - data->GetParameterValue(&value,gauss); - found=1; - break; - } - } - } - - if(found)*pvalue=value; - return found; } /*}}}*/ @@ -2346,11 +2108,4 @@ elementresult->ProcessUnits(this->parameters); } -} -/*}}}*/ -/*FUNCTION Tria::RheologyBbarx{{{1*/ -double Tria::RheologyBbarx(void){ - - return this->matice->GetBbar(); - } /*}}}*/ @@ -2598,6 +2353,10 @@ /*Get for Vx and Vy, the max of abs value: */ + #ifdef _HAVE_RESPONSES_ this->MaxAbsVx(&maxabsvx,false); this->MaxAbsVy(&maxabsvy,false); + #else + _error_("ISSM was not compiled with responses compiled in, exiting!"); + #endif /* Get node coordinates and dof list: */ @@ -2848,4 +2607,252 @@ } /*}}}*/ + +#ifdef _HAVE_RESPONSES_ +/*FUNCTION Tria::MassFlux {{{1*/ +double Tria::MassFlux( double* segment,bool process_units){ + + const int numdofs=2; + + int i; + double mass_flux=0; + double xyz_list[NUMVERTICES][3]; + double normal[2]; + double length,rho_ice; + double x1,y1,x2,y2,h1,h2; + double vx1,vx2,vy1,vy2; + GaussTria* gauss_1=NULL; + GaussTria* gauss_2=NULL; + + /*Get material parameters :*/ + rho_ice=matpar->GetRhoIce(); + + /*First off, check that this segment belongs to this element: */ + if ((int)*(segment+4)!=this->id)_error_("%s%i%s%i","error message: segment with id ",(int)*(segment+4)," does not belong to element with id:",this->id); + + /*Recover segment node locations: */ + x1=*(segment+0); y1=*(segment+1); x2=*(segment+2); y2=*(segment+3); + + /*Get xyz list: */ + GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES); + + /*get area coordinates of 0 and 1 locations: */ + gauss_1=new GaussTria(); + gauss_1->GaussFromCoords(x1,y1,&xyz_list[0][0]); + gauss_2=new GaussTria(); + gauss_2->GaussFromCoords(x2,y2,&xyz_list[0][0]); + + normal[0]=cos(atan2(x1-x2,y2-y1)); + normal[1]=sin(atan2(x1-x2,y2-y1)); + + length=sqrt(pow(x2-x1,2.0)+pow(y2-y1,2)); + + Input* thickness_input=inputs->GetInput(ThicknessEnum); _assert_(thickness_input); + Input* vx_input=inputs->GetInput(VxEnum); _assert_(vx_input); + Input* vy_input=inputs->GetInput(VyEnum); _assert_(vy_input); + + thickness_input->GetParameterValue(&h1, gauss_1); + thickness_input->GetParameterValue(&h2, gauss_2); + vx_input->GetParameterValue(&vx1,gauss_1); + vx_input->GetParameterValue(&vx2,gauss_2); + vy_input->GetParameterValue(&vy1,gauss_1); + vy_input->GetParameterValue(&vy2,gauss_2); + + mass_flux= rho_ice*length*( + (ONETHIRD*(h1-h2)*(vx1-vx2)+0.5*h2*(vx1-vx2)+0.5*(h1-h2)*vx2+h2*vx2)*normal[0]+ + (ONETHIRD*(h1-h2)*(vy1-vy2)+0.5*h2*(vy1-vy2)+0.5*(h1-h2)*vy2+h2*vy2)*normal[1] + ); + + /*Process units: */ + mass_flux=UnitConversion(mass_flux,IuToExtEnum,MassFluxEnum); + + /*clean up and return:*/ + delete gauss_1; + delete gauss_2; + return mass_flux; +} +/*}}}*/ +/*FUNCTION Tria::MaxAbsVx{{{1*/ +void Tria::MaxAbsVx(double* pmaxabsvx, bool process_units){ + + /*Get maximum:*/ + double maxabsvx=this->inputs->MaxAbs(VxEnum); + + /*process units if requested: */ + if(process_units) maxabsvx=UnitConversion(maxabsvx,IuToExtEnum,VxEnum); + + /*Assign output pointers:*/ + *pmaxabsvx=maxabsvx; +} +/*}}}*/ +/*FUNCTION Tria::MaxAbsVy{{{1*/ +void Tria::MaxAbsVy(double* pmaxabsvy, bool process_units){ + + /*Get maximum:*/ + double maxabsvy=this->inputs->MaxAbs(VyEnum); + + /*process units if requested: */ + if(process_units) maxabsvy=UnitConversion(maxabsvy,IuToExtEnum,VyEnum); + + /*Assign output pointers:*/ + *pmaxabsvy=maxabsvy; +} +/*}}}*/ +/*FUNCTION Tria::MaxAbsVz{{{1*/ +void Tria::MaxAbsVz(double* pmaxabsvz, bool process_units){ + + /*Get maximum:*/ + double maxabsvz=this->inputs->MaxAbs(VzEnum); + + /*process units if requested: */ + if(process_units) maxabsvz=UnitConversion(maxabsvz,IuToExtEnum,VyEnum); + + /*Assign output pointers:*/ + *pmaxabsvz=maxabsvz; +} +/*}}}*/ +/*FUNCTION Tria::MaxVel{{{1*/ +void Tria::MaxVel(double* pmaxvel, bool process_units){ + + /*Get maximum:*/ + double maxvel=this->inputs->Max(VelEnum); + + /*process units if requested: */ + if(process_units) maxvel=UnitConversion(maxvel,IuToExtEnum,VelEnum); + + /*Assign output pointers:*/ + *pmaxvel=maxvel; +} +/*}}}*/ +/*FUNCTION Tria::MaxVx{{{1*/ +void Tria::MaxVx(double* pmaxvx, bool process_units){ + + /*Get maximum:*/ + double maxvx=this->inputs->Max(VxEnum); + + /*process units if requested: */ + if(process_units) maxvx=UnitConversion(maxvx,IuToExtEnum,VxEnum); + + /*Assign output pointers:*/ + *pmaxvx=maxvx; +} +/*}}}*/ +/*FUNCTION Tria::MaxVy{{{1*/ +void Tria::MaxVy(double* pmaxvy, bool process_units){ + + /*Get maximum:*/ + double maxvy=this->inputs->Max(VyEnum); + + /*process units if requested: */ + if(process_units) maxvy=UnitConversion(maxvy,IuToExtEnum,VyEnum); + + /*Assign output pointers:*/ + *pmaxvy=maxvy; + +} +/*}}}*/ +/*FUNCTION Tria::MaxVz{{{1*/ +void Tria::MaxVz(double* pmaxvz, bool process_units){ + + /*Get maximum:*/ + double maxvz=this->inputs->Max(VzEnum); + + /*process units if requested: */ + if(process_units) maxvz=UnitConversion(maxvz,IuToExtEnum,VzEnum); + + /*Assign output pointers:*/ + *pmaxvz=maxvz; +} +/*}}}*/ +/*FUNCTION Tria::MinVel{{{1*/ +void Tria::MinVel(double* pminvel, bool process_units){ + + /*Get minimum:*/ + double minvel=this->inputs->Min(VelEnum); + + /*process units if requested: */ + if(process_units) minvel=UnitConversion(minvel,IuToExtEnum,VelEnum); + + /*Assign output pointers:*/ + *pminvel=minvel; +} +/*}}}*/ +/*FUNCTION Tria::MinVx{{{1*/ +void Tria::MinVx(double* pminvx, bool process_units){ + + /*Get minimum:*/ + double minvx=this->inputs->Min(VxEnum); + + /*process units if requested: */ + if(process_units) minvx=UnitConversion(minvx,IuToExtEnum,VxEnum); + + /*Assign output pointers:*/ + *pminvx=minvx; +} +/*}}}*/ +/*FUNCTION Tria::MinVy{{{1*/ +void Tria::MinVy(double* pminvy, bool process_units){ + + /*Get minimum:*/ + double minvy=this->inputs->Min(VyEnum); + + /*process units if requested: */ + if(process_units) minvy=UnitConversion(minvy,IuToExtEnum,VyEnum); + + /*Assign output pointers:*/ + *pminvy=minvy; +} +/*}}}*/ +/*FUNCTION Tria::MinVz{{{1*/ +void Tria::MinVz(double* pminvz, bool process_units){ + + /*Get minimum:*/ + double minvz=this->inputs->Min(VzEnum); + + /*process units if requested: */ + if(process_units) minvz=UnitConversion(minvz,IuToExtEnum,VzEnum); + + /*Assign output pointers:*/ + *pminvz=minvz; +} +/*}}}*/ +/*FUNCTION Tria::NodalValue {{{1*/ +int Tria::NodalValue(double* pvalue, int index, int natureofdataenum,bool process_units){ + + int i; + int found=0; + double value; + Input* data=NULL; + GaussTria *gauss = NULL; + + /*First, serarch the input: */ + data=inputs->GetInput(natureofdataenum); + + /*figure out if we have the vertex id: */ + found=0; + for(i=0;iGetVertexId()){ + /*Do we have natureofdataenum in our inputs? :*/ + if(data){ + /*ok, we are good. retrieve value of input at vertex :*/ + gauss=new GaussTria(); gauss->GaussVertex(i); + data->GetParameterValue(&value,gauss); + found=1; + break; + } + } + } + + if(found)*pvalue=value; + return found; +} +/*}}}*/ +/*FUNCTION Tria::RheologyBbarx{{{1*/ +double Tria::RheologyBbarx(void){ + + return this->matice->GetBbar(); + +} +/*}}}*/ +#endif #ifdef _HAVE_DIAGNOSTIC_ Index: /issm/trunk/src/c/objects/Elements/Tria.h =================================================================== --- /issm/trunk/src/c/objects/Elements/Tria.h (revision 9816) +++ /issm/trunk/src/c/objects/Elements/Tria.h (revision 9817) @@ -102,12 +102,5 @@ void DeleteResults(void); void MaterialUpdateFromTemperature(void){_error_("not implemented yet");}; - double MassFlux(double* segment,bool process_units); - void MaxAbsVx(double* pmaxabsvx, bool process_units); - void MaxAbsVy(double* pmaxabsvy, bool process_units); - void MaxAbsVz(double* pmaxabsvz, bool process_units); - void MaxVel(double* pmaxvel, bool process_units); - void MaxVx(double* pmaxvx, bool process_units); - void MaxVy(double* pmaxvy, bool process_units); - void MaxVz(double* pmaxvz, bool process_units); + void AgressiveMigration(); void SoftMigration(double* sheet_ungrounding); @@ -115,8 +108,5 @@ void PotentialSheetUngrounding(Vec potential_sheet_ungrounding); void MigrateGroundingLine(); - void MinVel(double* pminvel, bool process_units); - void MinVx(double* pminvx, bool process_units); - void MinVy(double* pminvy, bool process_units); - void MinVz(double* pminvz, bool process_units); + void RequestedOutput(int output_enum,int step,double time); void PatchFill(int* pcount, Patch* patch); @@ -130,6 +120,21 @@ double TimeAdapt(); int* GetHorizontalNeighboorSids(void); + + #ifdef _HAVE_RESPONSES_ + void MinVel(double* pminvel, bool process_units); + void MinVx(double* pminvx, bool process_units); + void MinVy(double* pminvy, bool process_units); + void MinVz(double* pminvz, bool process_units); double RheologyBbarx(void); int NodalValue(double* pvalue, int index, int natureofdataenum,bool process_units); + double MassFlux(double* segment,bool process_units); + void MaxAbsVx(double* pmaxabsvx, bool process_units); + void MaxAbsVy(double* pmaxabsvy, bool process_units); + void MaxAbsVz(double* pmaxabsvz, bool process_units); + void MaxVel(double* pmaxvel, bool process_units); + void MaxVx(double* pmaxvx, bool process_units); + void MaxVy(double* pmaxvy, bool process_units); + void MaxVz(double* pmaxvz, bool process_units); + #endif