Index: /u/astrid-r1b/morlighe/issmuci/trunk-jpl/../trunk-jpl/src/c/objects/Elements/Element.h =================================================================== --- /u/astrid-r1b/morlighe/issmuci/trunk-jpl/../trunk-jpl/src/c/objects/Elements/Element.h (revision 12470) +++ /u/astrid-r1b/morlighe/issmuci/trunk-jpl/../trunk-jpl/src/c/objects/Elements/Element.h (revision 12471) @@ -37,78 +37,78 @@ virtual int Sid()=0; virtual bool IsFloating()=0; virtual bool IsNodeOnShelf()=0; - virtual bool IsNodeOnShelfFromFlags(double* flags)=0; + virtual bool IsNodeOnShelfFromFlags(IssmDouble* flags)=0; virtual bool IsOnBed()=0; - virtual void GetInputListOnVertices(double* pvalue,int enumtype)=0; - virtual void GetInputListOnVertices(double* pvalue,int enumtype,double defaultvalue)=0; - virtual void GetInputValue(double* pvalue,Node* node,int enumtype)=0; + virtual void GetInputListOnVertices(IssmDouble* pvalue,int enumtype)=0; + virtual void GetInputListOnVertices(IssmDouble* pvalue,int enumtype,IssmDouble defaultvalue)=0; + virtual void GetInputValue(IssmDouble* pvalue,Node* node,int enumtype)=0; - virtual double SurfaceArea(void)=0; + virtual IssmDouble SurfaceArea(void)=0; virtual void InputDepthAverageAtBase(int enum_type,int average_enum_type,int object_enum)=0; virtual void ComputeBasalStress(Vector* sigma_b)=0; virtual void ComputeStrainRate(Vector* eps)=0; virtual void PatchSize(int* pnumrows, int* pnumvertices,int* pnumnodes)=0; virtual void PatchFill(int* pcount, Patch* patch)=0; - virtual void ListResultsInfo(int** results_enums,int** results_size,double** results_times,int** results_steps,int* num_results)=0; + virtual void ListResultsInfo(int** results_enums,int** results_size,IssmDouble** results_times,int** results_steps,int* num_results)=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 InputToResult(int enum_type,int step,IssmDouble 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 InputCreate(IssmDouble scalar,int name,int code)=0; + virtual void InputCreate(IssmDouble* 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 RequestedOutput(int output_enum,int step,IssmDouble time)=0; - virtual int NodalValue(double* pvalue, int index, int natureofdataenum,bool process_units)=0; - virtual void InputScale(int enum_type,double scale_factor)=0; + virtual int NodalValue(IssmDouble* pvalue, int index, int natureofdataenum,bool process_units)=0; + virtual void InputScale(int enum_type,IssmDouble scale_factor)=0; virtual void GetVectorFromInputs(Vector* vector, int name_enum)=0; virtual void GetVectorFromResults(Vector* vector,int id,int interp)=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(Vector* partition_contributions,Vector* partition_areas,double* vertex_response,double* qmu_part)=0; + virtual void InputArtificialNoise(int enum_type,IssmDouble min,IssmDouble max)=0; + virtual bool InputConvergence(IssmDouble* eps, int* enums,int num_enums,int* criterionenums,IssmDouble* criterionvalues,int num_criterionenums)=0; + virtual void AverageOntoPartition(Vector* partition_contributions,Vector* partition_areas,IssmDouble* vertex_response,IssmDouble* qmu_part)=0; virtual int* GetHorizontalNeighboorSids(void)=0; - virtual double TimeAdapt()=0; - virtual void MigrateGroundingLine(double* old_floating_ice,double* sheet_ungrounding)=0; + virtual IssmDouble TimeAdapt()=0; + virtual void MigrateGroundingLine(IssmDouble* old_floating_ice,IssmDouble* sheet_ungrounding)=0; virtual void PotentialSheetUngrounding(Vector* potential_sheet_ungrounding)=0; - virtual void PositiveDegreeDay(double* pdds,double* pds,double signorm)=0; - virtual int UpdatePotentialSheetUngrounding(double* potential_sheet_ungrounding,Vector* vec_nodes_on_iceshelf,double* nodes_on_iceshelf)=0; + virtual void PositiveDegreeDay(IssmDouble* pdds,IssmDouble* pds,IssmDouble signorm)=0; + virtual int UpdatePotentialSheetUngrounding(IssmDouble* potential_sheet_ungrounding,Vector* vec_nodes_on_iceshelf,IssmDouble* nodes_on_iceshelf)=0; virtual void ResetCoordinateSystem()=0; - virtual void SmearFunction(Vector* smearedvector,double (*WeightFunction)(double distance,double radius),double radius)=0; + virtual void SmearFunction(Vector* smearedvector,IssmDouble (*WeightFunction)(IssmDouble distance,IssmDouble radius),IssmDouble radius)=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 void ElementResponse(double* presponse,int response_enum,bool process_units)=0; - virtual double IceVolume(void)=0; + virtual void MinVel(IssmDouble* pminvel, bool process_units)=0; + virtual void MaxVel(IssmDouble* pmaxvel, bool process_units)=0; + virtual void MinVx(IssmDouble* pminvx, bool process_units)=0; + virtual void MaxVx(IssmDouble* pmaxvx, bool process_units)=0; + virtual void MaxAbsVx(IssmDouble* pmaxabsvx, bool process_units)=0; + virtual void MinVy(IssmDouble* pminvy, bool process_units)=0; + virtual void MaxVy(IssmDouble* pmaxvy, bool process_units)=0; + virtual void MaxAbsVy(IssmDouble* pmaxabsvy, bool process_units)=0; + virtual void MinVz(IssmDouble* pminvz, bool process_units)=0; + virtual void MaxVz(IssmDouble* pmaxvz, bool process_units)=0; + virtual void MaxAbsVz(IssmDouble* pmaxabsvz, bool process_units)=0; + virtual IssmDouble MassFlux(IssmDouble* segment,bool process_units)=0; + virtual void ElementResponse(IssmDouble* presponse,int response_enum,bool process_units)=0; + virtual IssmDouble IceVolume(void)=0; #endif #ifdef _HAVE_CONTROL_ virtual void Gradj(Vector* gradient,int control_type,int control_index)=0; - virtual double ThicknessAbsMisfit(bool process_units ,int weight_index)=0; - virtual double SurfaceAbsVelMisfit(bool process_units ,int weight_index)=0; - virtual double SurfaceRelVelMisfit(bool process_units ,int weight_index)=0; - virtual double SurfaceLogVelMisfit(bool process_units ,int weight_index)=0; - virtual double SurfaceLogVxVyMisfit(bool process_units,int weight_index)=0; - virtual double SurfaceAverageVelMisfit(bool process_units,int weight_index)=0; - virtual double ThicknessAbsGradient(bool process_units,int weight_index)=0; - virtual double RheologyBbarAbsGradient(bool process_units,int weight_index)=0; - virtual double DragCoefficientAbsGradient(bool process_units,int weight_index)=0; + virtual IssmDouble ThicknessAbsMisfit(bool process_units ,int weight_index)=0; + virtual IssmDouble SurfaceAbsVelMisfit(bool process_units ,int weight_index)=0; + virtual IssmDouble SurfaceRelVelMisfit(bool process_units ,int weight_index)=0; + virtual IssmDouble SurfaceLogVelMisfit(bool process_units ,int weight_index)=0; + virtual IssmDouble SurfaceLogVxVyMisfit(bool process_units,int weight_index)=0; + virtual IssmDouble SurfaceAverageVelMisfit(bool process_units,int weight_index)=0; + virtual IssmDouble ThicknessAbsGradient(bool process_units,int weight_index)=0; + virtual IssmDouble RheologyBbarAbsGradient(bool process_units,int weight_index)=0; + virtual IssmDouble DragCoefficientAbsGradient(bool process_units,int weight_index)=0; virtual void ControlInputGetGradient(Vector* gradient,int enum_type,int control_index)=0; - virtual void ControlInputSetGradient(double* gradient,int enum_type,int control_index)=0; - virtual void ControlInputScaleGradient(int enum_type, double scale)=0; + virtual void ControlInputSetGradient(IssmDouble* gradient,int enum_type,int control_index)=0; + virtual void ControlInputScaleGradient(int enum_type, IssmDouble scale)=0; virtual void GetVectorFromControlInputs(Vector* gradient,int control_enum,int control_index,const char* data)=0; - virtual void SetControlInputsFromVector(double* vector,int control_enum,int control_index)=0; - virtual void InputControlUpdate(double scalar,bool save_parameter)=0; + virtual void SetControlInputsFromVector(IssmDouble* vector,int control_enum,int control_index)=0; + virtual void InputControlUpdate(IssmDouble scalar,bool save_parameter)=0; #endif }; #endif Index: /u/astrid-r1b/morlighe/issmuci/trunk-jpl/../trunk-jpl/src/c/objects/Elements/Tria.cpp =================================================================== --- /u/astrid-r1b/morlighe/issmuci/trunk-jpl/../trunk-jpl/src/c/objects/Elements/Tria.cpp (revision 12470) +++ /u/astrid-r1b/morlighe/issmuci/trunk-jpl/../trunk-jpl/src/c/objects/Elements/Tria.cpp (revision 12471) @@ -124,16 +124,16 @@ /*Other*/ /*FUNCTION Tria::AverageOntoPartition {{{*/ -void Tria::AverageOntoPartition(Vector* partition_contributions,Vector* partition_areas,double* vertex_response,double* qmu_part){ +void Tria::AverageOntoPartition(Vector* partition_contributions,Vector* partition_areas,IssmDouble* vertex_response,IssmDouble* qmu_part){ bool already=false; int i,j; int partition[NUMVERTICES]; int offsetsid[NUMVERTICES]; int offsetdof[NUMVERTICES]; - double area; - double mean; - double values[3]; + IssmDouble area; + IssmDouble mean; + IssmDouble values[3]; /*First, get the area: */ area=this->GetArea(); @@ -231,10 +231,10 @@ /*Intermediaries */ int i,j,ig; - double heatcapacity,latentheat; - double Jdet,D_scalar; - double xyz_list[NUMVERTICES][3]; - double L[3]; + IssmDouble heatcapacity,latentheat; + IssmDouble Jdet,D_scalar; + IssmDouble xyz_list[NUMVERTICES][3]; + IssmDouble L[3]; GaussTria *gauss=NULL; /*Initialize Element matrix*/ @@ -290,18 +290,18 @@ /*Intermediaries */ int stabilization; int i,j,ig,dim; - double Jdettria,DL_scalar,dt,h; - double vel,vx,vy,dvxdx,dvydy; - double dvx[2],dvy[2]; - double v_gauss[2]={0.0}; - double xyz_list[NUMVERTICES][3]; - double L[NUMVERTICES]; - double B[2][NUMVERTICES]; - double Bprime[2][NUMVERTICES]; - double K[2][2] ={0.0}; - double KDL[2][2] ={0.0}; - double DL[2][2] ={0.0}; - double DLprime[2][2] ={0.0}; + IssmDouble Jdettria,DL_scalar,dt,h; + IssmDouble vel,vx,vy,dvxdx,dvydy; + IssmDouble dvx[2],dvy[2]; + IssmDouble v_gauss[2]={0.0}; + IssmDouble xyz_list[NUMVERTICES][3]; + IssmDouble L[NUMVERTICES]; + IssmDouble B[2][NUMVERTICES]; + IssmDouble Bprime[2][NUMVERTICES]; + IssmDouble K[2][2] ={0.0}; + IssmDouble KDL[2][2] ={0.0}; + IssmDouble DL[2][2] ={0.0}; + IssmDouble DLprime[2][2] ={0.0}; GaussTria *gauss=NULL; /*Initialize Element matrix*/ @@ -409,14 +409,14 @@ /*Intermediaries */ int i,j,ig,dim; - double xyz_list[NUMVERTICES][3]; - double Jdettria,dt,vx,vy; - double L[NUMVERTICES]; - double B[2][NUMVERTICES]; - double Bprime[2][NUMVERTICES]; - double DL[2][2]={0.0}; - double DLprime[2][2]={0.0}; - double DL_scalar; + IssmDouble xyz_list[NUMVERTICES][3]; + IssmDouble Jdettria,dt,vx,vy; + IssmDouble L[NUMVERTICES]; + IssmDouble B[2][NUMVERTICES]; + IssmDouble Bprime[2][NUMVERTICES]; + IssmDouble DL[2][2]={0.0}; + IssmDouble DLprime[2][2]={0.0}; + IssmDouble DL_scalar; GaussTria *gauss=NULL; /*Initialize Element matrix*/ @@ -484,9 +484,9 @@ /* Intermediaries */ int i,j,ig; - double DL_scalar,Jdet; - double xyz_list[NUMVERTICES][3]; - double L[1][3]; + IssmDouble DL_scalar,Jdet; + IssmDouble xyz_list[NUMVERTICES][3]; + IssmDouble L[1][3]; GaussTria *gauss = NULL; /*Initialize Element matrix*/ @@ -598,10 +598,10 @@ /*Intermediaries */ int i,j,ig; - double Jdettria,dt; - double surface_mass_balance_g,basal_melting_g,basal_melting_correction_g,thickness_g; - double xyz_list[NUMVERTICES][3]; - double L[NUMVERTICES]; + IssmDouble Jdettria,dt; + IssmDouble surface_mass_balance_g,basal_melting_g,basal_melting_correction_g,thickness_g; + IssmDouble xyz_list[NUMVERTICES][3]; + IssmDouble L[NUMVERTICES]; GaussTria* gauss=NULL; /*Initialize Element vector*/ @@ -649,10 +649,10 @@ /*Intermediaries */ int i,j,ig; - double Jdettria,dt; - double surface_mass_balance_g,basal_melting_g,thickness_g; - double xyz_list[NUMVERTICES][3]; - double L[NUMVERTICES]; + IssmDouble Jdettria,dt; + IssmDouble surface_mass_balance_g,basal_melting_g,thickness_g; + IssmDouble xyz_list[NUMVERTICES][3]; + IssmDouble L[NUMVERTICES]; GaussTria* gauss=NULL; /*Initialize Element vector*/ @@ -695,10 +695,10 @@ /*Intermediaries */ int i,j,ig; int analysis_type; - double Jdet; - double xyz_list[NUMVERTICES][3]; - double slope[2]; - double basis[3]; + IssmDouble Jdet; + IssmDouble xyz_list[NUMVERTICES][3]; + IssmDouble slope[2]; + IssmDouble basis[3]; GaussTria* gauss=NULL; /*Initialize Element vector*/ @@ -786,15 +786,15 @@ void Tria::ComputeStressTensor(){ int iv; - double xyz_list[NUMVERTICES][3]; - double pressure,viscosity; - double epsilon[3]; /* epsilon=[exx,eyy,exy];*/ - double sigma_xx[NUMVERTICES]; - double sigma_yy[NUMVERTICES]; - double sigma_zz[NUMVERTICES]={0,0,0}; - double sigma_xy[NUMVERTICES]; - double sigma_xz[NUMVERTICES]={0,0,0}; - double sigma_yz[NUMVERTICES]={0,0,0}; + IssmDouble xyz_list[NUMVERTICES][3]; + IssmDouble pressure,viscosity; + IssmDouble epsilon[3]; /* epsilon=[exx,eyy,exy];*/ + IssmDouble sigma_xx[NUMVERTICES]; + IssmDouble sigma_yy[NUMVERTICES]; + IssmDouble sigma_zz[NUMVERTICES]={0,0,0}; + IssmDouble sigma_xy[NUMVERTICES]; + IssmDouble sigma_xz[NUMVERTICES]={0,0,0}; + IssmDouble sigma_yz[NUMVERTICES]={0,0,0}; GaussTria* gauss=NULL; /* Get node coordinates and dof list: */ @@ -947,11 +947,11 @@ } /*}}}*/ /*FUNCTION Tria::GetArea {{{*/ -double Tria::GetArea(void){ +IssmDouble Tria::GetArea(void){ - double area=0; - double xyz_list[NUMVERTICES][3]; - double x1,y1,x2,y2,x3,y3; + IssmDouble area=0; + IssmDouble xyz_list[NUMVERTICES][3]; + IssmDouble x1,y1,x2,y2,x3,y3; /*Get xyz list: */ GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES); @@ -1021,11 +1021,11 @@ _error_("Node provided not found among element nodes"); } /*}}}*/ -/*FUNCTION Tria::GetInputListOnVertices(double* pvalue,int enumtype) {{{*/ -void Tria::GetInputListOnVertices(double* pvalue,int enumtype){ +/*FUNCTION Tria::GetInputListOnVertices(IssmDouble* pvalue,int enumtype) {{{*/ +void Tria::GetInputListOnVertices(IssmDouble* pvalue,int enumtype){ /*Intermediaries*/ - double value[NUMVERTICES]; + IssmDouble value[NUMVERTICES]; GaussTria *gauss = NULL; /*Recover input*/ @@ -1046,10 +1046,10 @@ delete gauss; } /*}}}*/ -/*FUNCTION Tria::GetInputListOnVertices(double* pvalue,int enumtype,double defaultvalue) {{{*/ -void Tria::GetInputListOnVertices(double* pvalue,int enumtype,double defaultvalue){ +/*FUNCTION Tria::GetInputListOnVertices(IssmDouble* pvalue,int enumtype,IssmDouble defaultvalue) {{{*/ +void Tria::GetInputListOnVertices(IssmDouble* pvalue,int enumtype,IssmDouble defaultvalue){ - double value[NUMVERTICES]; + IssmDouble value[NUMVERTICES]; GaussTria *gauss = NULL; Input *input = inputs->GetInput(enumtype); @@ -1072,10 +1072,10 @@ delete gauss; } /*}}}*/ -/*FUNCTION Tria::GetInputListOnVertices(double* pvalue,int enumtype,double defaultvalue,int index) TO BE REMOVED{{{*/ -void Tria::GetInputListOnVertices(double* pvalue,int enumtype,double defaultvalue,int index){ +/*FUNCTION Tria::GetInputListOnVertices(IssmDouble* pvalue,int enumtype,IssmDouble defaultvalue,int index) TO BE REMOVED{{{*/ +void Tria::GetInputListOnVertices(IssmDouble* pvalue,int enumtype,IssmDouble defaultvalue,int index){ - double value[NUMVERTICES]; + IssmDouble value[NUMVERTICES]; GaussTria *gauss = NULL; Input *input = inputs->GetInput(enumtype); @@ -1098,8 +1098,8 @@ delete gauss; } /*}}}*/ -/*FUNCTION Tria::GetInputValue(double* pvalue,Node* node,int enumtype) {{{*/ -void Tria::GetInputValue(double* pvalue,Node* node,int enumtype){ +/*FUNCTION Tria::GetInputValue(IssmDouble* pvalue,Node* node,int enumtype) {{{*/ +void Tria::GetInputValue(IssmDouble* pvalue,Node* node,int enumtype){ Input* input=inputs->GetInput(enumtype); if(!input) _error_("No input of type %s found in tria",EnumToStringx(enumtype)); @@ -1149,14 +1149,14 @@ } /*}}}*/ -/*FUNCTION Tria::GetStrainRate2d(double* epsilon,double* xyz_list, GaussTria* gauss, Input* vx_input, Input* vy_input){{{*/ -void Tria::GetStrainRate2d(double* epsilon,double* xyz_list, GaussTria* gauss, Input* vx_input, Input* vy_input){ +/*FUNCTION Tria::GetStrainRate2d(IssmDouble* epsilon,IssmDouble* xyz_list, GaussTria* gauss, Input* vx_input, Input* vy_input){{{*/ +void Tria::GetStrainRate2d(IssmDouble* epsilon,IssmDouble* xyz_list, GaussTria* gauss, Input* vx_input, Input* vy_input){ /*Compute the 2d Strain Rate (3 components): * epsilon=[exx eyy exy] */ int i; - double epsilonvx[3]; - double epsilonvy[3]; + IssmDouble epsilonvx[3]; + IssmDouble epsilonvy[3]; /*Check that both inputs have been found*/ if (!vx_input || !vy_input){ @@ -1225,7 +1225,7 @@ } /*}}}*/ /*FUNCTION Tria::InputArtificialNoise{{{*/ -void Tria::InputArtificialNoise(int enum_type,double min,double max){ +void Tria::InputArtificialNoise(int enum_type,IssmDouble min,IssmDouble max){ Input* input=NULL; @@ -1238,7 +1238,7 @@ } /*}}}*/ /*FUNCTION Tria::InputConvergence{{{*/ -bool Tria::InputConvergence(double* eps, int* enums,int num_enums,int* criterionenums,double* criterionvalues,int num_criterionenums){ +bool Tria::InputConvergence(IssmDouble* eps, int* enums,int num_enums,int* criterionenums,IssmDouble* criterionvalues,int num_criterionenums){ bool converged=true; int i; @@ -1305,7 +1305,7 @@ } /*}}}*/ /*FUNCTION Tria::InputScale{{{*/ -void Tria::InputScale(int enum_type,double scale_factor){ +void Tria::InputScale(int enum_type,IssmDouble scale_factor){ Input* input=NULL; @@ -1318,7 +1318,7 @@ } /*}}}*/ /*FUNCTION Tria::InputToResult{{{*/ -void Tria::InputToResult(int enum_type,int step,double time){ +void Tria::InputToResult(int enum_type,int step,IssmDouble time){ int i; Input *input = NULL; @@ -1349,8 +1349,8 @@ this->inputs->AddInput(new IntInput(name,constant)); } /*}}}*/ -/*FUNCTION Tria::InputUpdateFromConstant(double value, int name);{{{*/ -void Tria::InputUpdateFromConstant(double constant, int name){ +/*FUNCTION Tria::InputUpdateFromConstant(IssmDouble value, int name);{{{*/ +void Tria::InputUpdateFromConstant(IssmDouble constant, int name){ /*Check that name is an element input*/ if (!IsInput(name)) return; @@ -1373,12 +1373,12 @@ /*Intermediaries*/ int i,j; int tria_vertex_ids[3]; - double nodeinputs[3]; - double cmmininputs[3]; - double cmmaxinputs[3]; + IssmDouble nodeinputs[3]; + IssmDouble cmmininputs[3]; + IssmDouble cmmaxinputs[3]; bool control_analysis=false; int num_control_type; - double yts; + IssmDouble yts; int num_cm_responses; /*Get parameters: */ @@ -1454,7 +1454,7 @@ } /*}}}*/ /*FUNCTION Tria::InputUpdateFromSolution {{{*/ -void Tria::InputUpdateFromSolution(double* solution){ +void Tria::InputUpdateFromSolution(IssmDouble* solution){ /*retrive parameters: */ int analysis_type; @@ -1509,12 +1509,12 @@ } /*}}}*/ /*FUNCTION Tria::InputUpdateFromSolutionOneDof{{{*/ -void Tria::InputUpdateFromSolutionOneDof(double* solution,int enum_type){ +void Tria::InputUpdateFromSolutionOneDof(IssmDouble* solution,int enum_type){ const int numdof = NDOF1*NUMVERTICES; int* doflist=NULL; - double values[numdof]; + IssmDouble values[numdof]; /*Get dof list: */ GetDofList(&doflist,NoneApproximationEnum,GsetEnum); @@ -1533,20 +1533,20 @@ } /*}}}*/ /*FUNCTION Tria::InputUpdateFromSolutionPrognostic{{{*/ -void Tria::InputUpdateFromSolutionPrognostic(double* solution){ +void Tria::InputUpdateFromSolutionPrognostic(IssmDouble* solution){ /*Intermediaries*/ const int numdof = NDOF1*NUMVERTICES; int i,hydroadjustment; int* doflist=NULL; - double rho_ice,rho_water,minthickness; - double newthickness[numdof]; - double newbed[numdof]; - double newsurface[numdof]; - double oldbed[NUMVERTICES]; - double oldsurface[NUMVERTICES]; - double oldthickness[NUMVERTICES]; + IssmDouble rho_ice,rho_water,minthickness; + IssmDouble newthickness[numdof]; + IssmDouble newbed[numdof]; + IssmDouble newsurface[numdof]; + IssmDouble oldbed[NUMVERTICES]; + IssmDouble oldsurface[NUMVERTICES]; + IssmDouble oldthickness[NUMVERTICES]; /*Get dof list: */ GetDofList(&doflist,NoneApproximationEnum,GsetEnum); @@ -1599,8 +1599,8 @@ xDelete(doflist); } /*}}}*/ -/*FUNCTION Tria::InputUpdateFromVector(double* vector, int name, int type);{{{*/ -void Tria::InputUpdateFromVector(double* vector, int name, int type){ +/*FUNCTION Tria::InputUpdateFromVector(IssmDouble* vector, int name, int type);{{{*/ +void Tria::InputUpdateFromVector(IssmDouble* vector, int name, int type){ /*Check that name is an element input*/ if (!IsInput(name)) return; @@ -1610,7 +1610,7 @@ case VertexEnum: /*New TriaP1Input*/ - double values[3]; + IssmDouble values[3]; /*Get values on the 3 vertices*/ for (int i=0;i<3;i++){ @@ -1641,8 +1641,8 @@ _error_(" not supported yet!"); } /*}}}*/ -/*FUNCTION Tria::InputCreate(double scalar,int enum,int code);{{{*/ -void Tria::InputCreate(double scalar,int name,int code){ +/*FUNCTION Tria::InputCreate(IssmDouble scalar,int enum,int code);{{{*/ +void Tria::InputCreate(IssmDouble scalar,int name,int code){ /*Check that name is an element input*/ if (!IsInput(name)) return; @@ -1653,26 +1653,26 @@ else if ((code==6) || (code==2)){ //integer this->inputs->AddInput(new IntInput(name,(int)scalar)); } - else if ((code==7) || (code==3)){ //double + else if ((code==7) || (code==3)){ //IssmDouble this->inputs->AddInput(new DoubleInput(name,(int)scalar)); } else _error_("%s%i"," could not recognize nature of vector from code ",code); } /*}}}*/ -/*FUNCTION Tria::InputCreate(double* vector,int index,IoModel* iomodel,int M,int N,int vector_type,int vector_enum,int code){{{*/ -void Tria::InputCreate(double* vector, int index,IoModel* iomodel,int M,int N,int vector_type,int vector_enum,int code){//index into elements +/*FUNCTION Tria::InputCreate(IssmDouble* vector,int index,IoModel* iomodel,int M,int N,int vector_type,int vector_enum,int code){{{*/ +void Tria::InputCreate(IssmDouble* vector, int index,IoModel* iomodel,int M,int N,int vector_type,int vector_enum,int code){//index into elements /*Intermediaries*/ int i,j,t; int tria_vertex_ids[3]; int row; - double nodeinputs[3]; - double time; + IssmDouble nodeinputs[3]; + IssmDouble time; TransientInput* transientinput=NULL; int numberofvertices; int numberofelements; - double yts; + IssmDouble yts; /*Fetch parameters: */ @@ -1692,7 +1692,7 @@ if(M==numberofvertices){ /*create input values: */ - for(i=0;i<3;i++)nodeinputs[i]=(double)vector[tria_vertex_ids[i]-1]; + for(i=0;i<3;i++)nodeinputs[i]=(IssmDouble)vector[tria_vertex_ids[i]-1]; /*process units: */ UnitConversion(&nodeinputs[0], 3 ,ExtToIuEnum, vector_enum); @@ -1707,14 +1707,14 @@ /*create input values: */ for(i=0;i<3;i++){ row=tria_vertex_ids[i]-1; - nodeinputs[i]=(double)vector[N*row+t]; + nodeinputs[i]=(IssmDouble)vector[N*row+t]; } /*process units: */ UnitConversion(&nodeinputs[0], 3 ,ExtToIuEnum, vector_enum); /*time? :*/ - time=(double)vector[(M-1)*N+t]*yts; + time=(IssmDouble)vector[(M-1)*N+t]*yts; if(t==0) transientinput=new TransientInput(vector_enum); transientinput->AddTimeInput(new TriaP1Input(vector_enum,nodeinputs),time); @@ -1735,8 +1735,8 @@ else if (code==6){ //integer this->inputs->AddInput(new IntInput(vector_enum,(int)vector[index])); } - else if (code==7){ //double - this->inputs->AddInput(new DoubleInput(vector_enum,(double)vector[index])); + else if (code==7){ //IssmDouble + this->inputs->AddInput(new DoubleInput(vector_enum,(IssmDouble)vector[index])); } else _error_("%s%i"," could not recognize nature of vector from code ",code); } @@ -1816,7 +1816,7 @@ } /*}}}*/ /*FUNCTION Tria::IsNodeOnShelfFromFlags {{{*/ -bool Tria::IsNodeOnShelfFromFlags(double* flags){ +bool Tria::IsNodeOnShelfFromFlags(IssmDouble* flags){ int i; bool shelf=false; @@ -1839,14 +1839,14 @@ } /*}}}*/ /*FUNCTION Tria::ListResultsInfo{{{*/ -void Tria::ListResultsInfo(int** in_resultsenums,int** in_resultssizes,double** in_resultstimes,int** in_resultssteps,int* in_num_results){ +void Tria::ListResultsInfo(int** in_resultsenums,int** in_resultssizes,IssmDouble** in_resultstimes,int** in_resultssteps,int* in_num_results){ /*Intermediaries*/ int i; int numberofresults = 0; int *resultsenums = NULL; int *resultssizes = NULL; - double *resultstimes = NULL; + IssmDouble *resultstimes = NULL; int *resultssteps = NULL; /*Checks*/ @@ -1863,7 +1863,7 @@ /*Allocate output*/ resultsenums=xNew(numberofresults); resultssizes=xNew(numberofresults); - resultstimes=xNew(numberofresults); + resultstimes=xNew(numberofresults); resultssteps=xNew(numberofresults); /*populate enums*/ @@ -1890,14 +1890,14 @@ }/*}}}*/ /*FUNCTION Tria::MigrateGroundingLine{{{*/ -void Tria::MigrateGroundingLine(double* old_floating_ice,double* sheet_ungrounding){ +void Tria::MigrateGroundingLine(IssmDouble* old_floating_ice,IssmDouble* sheet_ungrounding){ int i,migration_style,unground; bool elementonshelf = false; - double bed_hydro,yts,gl_melting_rate; - double rho_water,rho_ice,density; - double melting[NUMVERTICES]; - double h[NUMVERTICES],s[NUMVERTICES],b[NUMVERTICES],ba[NUMVERTICES]; + IssmDouble bed_hydro,yts,gl_melting_rate; + IssmDouble rho_water,rho_ice,density; + IssmDouble melting[NUMVERTICES]; + IssmDouble h[NUMVERTICES],s[NUMVERTICES],b[NUMVERTICES],ba[NUMVERTICES]; /*Recover info at the vertices: */ parameters->FindParam(&migration_style,GroundinglineMigrationEnum); @@ -1972,11 +1972,11 @@ } /*}}}*/ /*FUNCTION Tria::NodalValue {{{*/ -int Tria::NodalValue(double* pvalue, int index, int natureofdataenum,bool process_units){ +int Tria::NodalValue(IssmDouble* pvalue, int index, int natureofdataenum,bool process_units){ int i; int found=0; - double value; + IssmDouble value; Input* data=NULL; GaussTria *gauss = NULL; @@ -2057,9 +2057,9 @@ void Tria::PotentialSheetUngrounding(Vector* potential_sheet_ungrounding){ int i; - double h[NUMVERTICES],ba[NUMVERTICES]; - double bed_hydro; - double rho_water,rho_ice,density; + IssmDouble h[NUMVERTICES],ba[NUMVERTICES]; + IssmDouble bed_hydro; + IssmDouble rho_water,rho_ice,density; bool elementonshelf = false; /*material parameters: */ @@ -2083,58 +2083,58 @@ } /*}}}*/ /*FUNCTION Tria::PositiveDegreeDay{{{*/ -void Tria::PositiveDegreeDay(double* pdds,double* pds,double signorm){ +void Tria::PositiveDegreeDay(IssmDouble* pdds,IssmDouble* pds,IssmDouble signorm){ int i,iqj,imonth; - double agd[NUMVERTICES]; // surface mass balance - double saccu[NUMVERTICES] = {0}; // yearly surface accumulation - double smelt[NUMVERTICES] = {0}; // yearly melt - double precrunoff[NUMVERTICES]; // yearly runoff - double prect; // total precipitation during 1 year taking into account des. ef. - double water; //water=rain + snowmelt - double runoff; //meltwater only, does not include rain - double sconv; //rhow_rain/rhoi / 12 months + IssmDouble agd[NUMVERTICES]; // surface mass balance + IssmDouble saccu[NUMVERTICES] = {0}; // yearly surface accumulation + IssmDouble smelt[NUMVERTICES] = {0}; // yearly melt + IssmDouble precrunoff[NUMVERTICES]; // yearly runoff + IssmDouble prect; // total precipitation during 1 year taking into account des. ef. + IssmDouble water; //water=rain + snowmelt + IssmDouble runoff; //meltwater only, does not include rain + IssmDouble sconv; //rhow_rain/rhoi / 12 months - double rho_water,rho_ice,density; - double lapser=6.5/1000, sealev=0; // lapse rate. degrees per meter. 7.5 lev's 99 paper, 9 Marshall 99 paper - double desfac = 0.5; // desert elevation factor - double s0p[NUMVERTICES]={0}; // should be set to elevation from precip source - double s0t[NUMVERTICES]={0}; // should be set to elevation from temperature source - double st; // elevation between altitude of the temp record and current altitude - double sp; // elevation between altitude of the prec record and current altitude + IssmDouble rho_water,rho_ice,density; + IssmDouble lapser=6.5/1000, sealev=0; // lapse rate. degrees per meter. 7.5 lev's 99 paper, 9 Marshall 99 paper + IssmDouble desfac = 0.5; // desert elevation factor + IssmDouble s0p[NUMVERTICES]={0}; // should be set to elevation from precip source + IssmDouble s0t[NUMVERTICES]={0}; // should be set to elevation from temperature source + IssmDouble st; // elevation between altitude of the temp record and current altitude + IssmDouble sp; // elevation between altitude of the prec record and current altitude // PDD and PD constants and variables - double siglim; // sigma limit for the integration which is equal to 2.5 sigmanorm - double signormc = signorm - 0.5; // sigma of the temperature distribution for cloudy day - double siglimc, siglim0, siglim0c; - double PDup, pddsig, PDCUT = 2.0; // PDcut: rain/snow cutoff temperature (C) - double DT = 0.02; - double pddt, pd; // pd: snow/precip fraction, precipitation falling as snow + IssmDouble siglim; // sigma limit for the integration which is equal to 2.5 sigmanorm + IssmDouble signormc = signorm - 0.5; // sigma of the temperature distribution for cloudy day + IssmDouble siglimc, siglim0, siglim0c; + IssmDouble PDup, pddsig, PDCUT = 2.0; // PDcut: rain/snow cutoff temperature (C) + IssmDouble DT = 0.02; + IssmDouble pddt, pd; // pd: snow/precip fraction, precipitation falling as snow - double q, qmpt; // q is desert/elev. fact, hnpfac is huybrect fact, and pd is normal dist. - double qm[NUMVERTICES] = {0}; // snow part of the precipitation - double qmt[NUMVERTICES] = {0}; // precipitation without desertification effect adjustment - double qmp[NUMVERTICES] = {0}; // desertification taken into account - double pdd[NUMVERTICES] = {0}; - double frzndd[NUMVERTICES] = {0}; + IssmDouble q, qmpt; // q is desert/elev. fact, hnpfac is huybrect fact, and pd is normal dist. + IssmDouble qm[NUMVERTICES] = {0}; // snow part of the precipitation + IssmDouble qmt[NUMVERTICES] = {0}; // precipitation without desertification effect adjustment + IssmDouble qmp[NUMVERTICES] = {0}; // desertification taken into account + IssmDouble pdd[NUMVERTICES] = {0}; + IssmDouble frzndd[NUMVERTICES] = {0}; - double tstar; // monthly mean surface temp - double Tsum[NUMVERTICES]= {0}; // average summer (JJA) temperature - double Tsurf[NUMVERTICES] = {0}; // average annual temperature + IssmDouble tstar; // monthly mean surface temp + IssmDouble Tsum[NUMVERTICES]= {0}; // average summer (JJA) temperature + IssmDouble Tsurf[NUMVERTICES] = {0}; // average annual temperature - double h[NUMVERTICES],s[NUMVERTICES],ttmp[NUMVERTICES],prectmp[NUMVERTICES]; // ,b[NUMVERTICES] - double t[NUMVERTICES+1][12],prec[NUMVERTICES+1][12]; - double deltm=1/12; + IssmDouble h[NUMVERTICES],s[NUMVERTICES],ttmp[NUMVERTICES],prectmp[NUMVERTICES]; // ,b[NUMVERTICES] + IssmDouble t[NUMVERTICES+1][12],prec[NUMVERTICES+1][12]; + IssmDouble deltm=1/12; int ismon[12]={12,1,2,3,4,5,6,7,8,9,10,11}; - double snwm; // snow that could have been melted in a year. - double snwmf; // ablation factor for snow per positive degree day. - double smf; // ablation factor for ice per pdd (Braithwaite 1995 from tarasov 2002). + IssmDouble snwm; // snow that could have been melted in a year. + IssmDouble snwmf; // ablation factor for snow per positive degree day. + IssmDouble smf; // ablation factor for ice per pdd (Braithwaite 1995 from tarasov 2002). - double dfrz=1.5, CovrLm=2009./3.35e+5, dCovrLm=dfrz*CovrLm; //m*J kg^-1 C^-1 /(J kg^-1)=m/C yr - double supice,supcap,diffndd; - double fsupT=0.5, fsupndd=0.5; // Tsurf mode factors for supice - double pddtj[NUMVERTICES], hmx2; + IssmDouble dfrz=1.5, CovrLm=2009./3.35e+5, dCovrLm=dfrz*CovrLm; //m*J kg^-1 C^-1 /(J kg^-1)=m/C yr + IssmDouble supice,supcap,diffndd; + IssmDouble fsupT=0.5, fsupndd=0.5; // Tsurf mode factors for supice + IssmDouble pddtj[NUMVERTICES], hmx2; /*Recover info at the vertices: */ GetInputListOnVertices(&h[0],ThicknessEnum); @@ -2145,7 +2145,7 @@ /*Recover monthly temperatures*/ Input* input=inputs->GetInput(SurfaceforcingsMonthlytemperaturesEnum); _assert_(input); GaussTria* gauss=new GaussTria(); - double time,yts; + IssmDouble time,yts; this->parameters->FindParam(&time,TimeEnum); this->parameters->FindParam(&yts,ConstantsYtsEnum); for(int month=0;month<12;month++){ @@ -2329,7 +2329,7 @@ } /*}}}*/ /*FUNCTION Tria::RequestedOutput{{{*/ -void Tria::RequestedOutput(int output_enum,int step,double time){ +void Tria::RequestedOutput(int output_enum,int step,IssmDouble time){ if(IsInput(output_enum)){ /*just transfer this input to results, and we are done: */ @@ -2363,7 +2363,7 @@ } /*}}}*/ /*FUNCTION Tria::SmearFunction {{{*/ -void Tria::SmearFunction(Vector* smearedvector,double (*WeightFunction)(double distance,double radius),double radius){ +void Tria::SmearFunction(Vector* smearedvector,IssmDouble (*WeightFunction)(IssmDouble distance,IssmDouble radius),IssmDouble radius){ _error_("not implemented yet"); } @@ -2385,13 +2385,13 @@ } /*}}}*/ /*FUNCTION Tria::SurfaceArea {{{*/ -double Tria::SurfaceArea(void){ +IssmDouble Tria::SurfaceArea(void){ int i; - double S; - double normal[3]; - double v13[3],v23[3]; - double xyz_list[NUMVERTICES][3]; + IssmDouble S; + IssmDouble normal[3]; + IssmDouble v13[3],v23[3]; + IssmDouble xyz_list[NUMVERTICES][3]; /*If on water, return 0: */ if(IsOnWater())return 0; @@ -2407,19 +2407,19 @@ normal[1]=v13[2]*v23[0]-v13[0]*v23[2]; normal[2]=v13[0]*v23[1]-v13[1]*v23[0]; - S = 0.5 * sqrt(pow(normal[0],(double)2)+pow(normal[1],(double)2)+pow(normal[2],(double)2)); + S = 0.5 * sqrt(pow(normal[0],(IssmDouble)2)+pow(normal[1],(IssmDouble)2)+pow(normal[2],(IssmDouble)2)); /*Return: */ return S; } /*}}}*/ /*FUNCTION Tria::SurfaceNormal{{{*/ -void Tria::SurfaceNormal(double* surface_normal, double xyz_list[3][3]){ +void Tria::SurfaceNormal(IssmDouble* surface_normal, IssmDouble xyz_list[3][3]){ int i; - double v13[3],v23[3]; - double normal[3]; - double normal_norm; + IssmDouble v13[3],v23[3]; + IssmDouble normal[3]; + IssmDouble normal_norm; for (i=0;i<3;i++){ v13[i]=xyz_list[0][i]-xyz_list[2][i]; @@ -2430,7 +2430,7 @@ normal[1]=v13[2]*v23[0]-v13[0]*v23[2]; normal[2]=v13[0]*v23[1]-v13[1]*v23[0]; - normal_norm=sqrt( pow(normal[0],(double)2)+pow(normal[1],(double)2)+pow(normal[2],(double)2) ); + normal_norm=sqrt( pow(normal[0],(IssmDouble)2)+pow(normal[1],(IssmDouble)2)+pow(normal[2],(IssmDouble)2) ); *(surface_normal)=normal[0]/normal_norm; *(surface_normal+1)=normal[1]/normal_norm; @@ -2438,16 +2438,16 @@ } /*}}}*/ /*FUNCTION Tria::TimeAdapt{{{*/ -double Tria::TimeAdapt(void){ +IssmDouble Tria::TimeAdapt(void){ /*intermediary: */ int i; - double C,dt; - double dx,dy; - double maxx,minx; - double maxy,miny; - double maxabsvx,maxabsvy; - double xyz_list[NUMVERTICES][3]; + IssmDouble C,dt; + IssmDouble dx,dy; + IssmDouble maxx,minx; + IssmDouble maxy,miny; + IssmDouble maxabsvx,maxabsvy; + IssmDouble xyz_list[NUMVERTICES][3]; /*get CFL coefficient:*/ this->parameters->FindParam(&C,TimesteppingCflCoefficientEnum); @@ -2491,8 +2491,8 @@ int tria_node_ids[3]; int tria_vertex_ids[3]; int tria_type; - double nodeinputs[3]; - double yts; + IssmDouble nodeinputs[3]; + IssmDouble yts; int progstabilization,balancestabilization; bool dakota_analysis; @@ -2584,7 +2584,7 @@ } /*}}}*/ /*FUNCTION Tria::UpdatePotentialSheetUngrounding{{{*/ -int Tria::UpdatePotentialSheetUngrounding(double* vertices_potentially_ungrounding,Vector* vec_nodes_on_iceshelf,double* nodes_on_iceshelf){ +int Tria::UpdatePotentialSheetUngrounding(IssmDouble* vertices_potentially_ungrounding,Vector* vec_nodes_on_iceshelf,IssmDouble* nodes_on_iceshelf){ int i; int nflipped=0; @@ -2606,11 +2606,11 @@ #ifdef _HAVE_RESPONSES_ /*FUNCTION Tria::IceVolume {{{*/ -double Tria::IceVolume(void){ +IssmDouble Tria::IceVolume(void){ /*The volume of a troncated prism is base * 1/3 sum(length of edges)*/ - double base,surface,bed; - double xyz_list[NUMVERTICES][3]; + IssmDouble base,surface,bed; + IssmDouble xyz_list[NUMVERTICES][3]; if(IsOnWater())return 0; @@ -2632,17 +2632,17 @@ } /*}}}*/ /*FUNCTION Tria::MassFlux {{{*/ -double Tria::MassFlux( double* segment,bool process_units){ +IssmDouble Tria::MassFlux( IssmDouble* segment,bool process_units){ const int numdofs=2; int i,dim; - 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; + IssmDouble mass_flux=0; + IssmDouble xyz_list[NUMVERTICES][3]; + IssmDouble normal[2]; + IssmDouble length,rho_ice; + IssmDouble x1,y1,x2,y2,h1,h2; + IssmDouble vx1,vx2,vy1,vy2; GaussTria* gauss_1=NULL; GaussTria* gauss_2=NULL; @@ -2704,10 +2704,10 @@ } /*}}}*/ /*FUNCTION Tria::MaxAbsVx{{{*/ -void Tria::MaxAbsVx(double* pmaxabsvx, bool process_units){ +void Tria::MaxAbsVx(IssmDouble* pmaxabsvx, bool process_units){ /*Get maximum:*/ - double maxabsvx=this->inputs->MaxAbs(VxEnum); + IssmDouble maxabsvx=this->inputs->MaxAbs(VxEnum); /*process units if requested: */ if(process_units) maxabsvx=UnitConversion(maxabsvx,IuToExtEnum,VxEnum); @@ -2717,10 +2717,10 @@ } /*}}}*/ /*FUNCTION Tria::MaxAbsVy{{{*/ -void Tria::MaxAbsVy(double* pmaxabsvy, bool process_units){ +void Tria::MaxAbsVy(IssmDouble* pmaxabsvy, bool process_units){ /*Get maximum:*/ - double maxabsvy=this->inputs->MaxAbs(VyEnum); + IssmDouble maxabsvy=this->inputs->MaxAbs(VyEnum); /*process units if requested: */ if(process_units) maxabsvy=UnitConversion(maxabsvy,IuToExtEnum,VyEnum); @@ -2730,10 +2730,10 @@ } /*}}}*/ /*FUNCTION Tria::MaxAbsVz{{{*/ -void Tria::MaxAbsVz(double* pmaxabsvz, bool process_units){ +void Tria::MaxAbsVz(IssmDouble* pmaxabsvz, bool process_units){ /*Get maximum:*/ - double maxabsvz=this->inputs->MaxAbs(VzEnum); + IssmDouble maxabsvz=this->inputs->MaxAbs(VzEnum); /*process units if requested: */ if(process_units) maxabsvz=UnitConversion(maxabsvz,IuToExtEnum,VyEnum); @@ -2743,10 +2743,10 @@ } /*}}}*/ /*FUNCTION Tria::MaxVel{{{*/ -void Tria::MaxVel(double* pmaxvel, bool process_units){ +void Tria::MaxVel(IssmDouble* pmaxvel, bool process_units){ /*Get maximum:*/ - double maxvel=this->inputs->Max(VelEnum); + IssmDouble maxvel=this->inputs->Max(VelEnum); /*process units if requested: */ if(process_units) maxvel=UnitConversion(maxvel,IuToExtEnum,VelEnum); @@ -2756,10 +2756,10 @@ } /*}}}*/ /*FUNCTION Tria::MaxVx{{{*/ -void Tria::MaxVx(double* pmaxvx, bool process_units){ +void Tria::MaxVx(IssmDouble* pmaxvx, bool process_units){ /*Get maximum:*/ - double maxvx=this->inputs->Max(VxEnum); + IssmDouble maxvx=this->inputs->Max(VxEnum); /*process units if requested: */ if(process_units) maxvx=UnitConversion(maxvx,IuToExtEnum,VxEnum); @@ -2769,10 +2769,10 @@ } /*}}}*/ /*FUNCTION Tria::MaxVy{{{*/ -void Tria::MaxVy(double* pmaxvy, bool process_units){ +void Tria::MaxVy(IssmDouble* pmaxvy, bool process_units){ /*Get maximum:*/ - double maxvy=this->inputs->Max(VyEnum); + IssmDouble maxvy=this->inputs->Max(VyEnum); /*process units if requested: */ if(process_units) maxvy=UnitConversion(maxvy,IuToExtEnum,VyEnum); @@ -2783,10 +2783,10 @@ } /*}}}*/ /*FUNCTION Tria::MaxVz{{{*/ -void Tria::MaxVz(double* pmaxvz, bool process_units){ +void Tria::MaxVz(IssmDouble* pmaxvz, bool process_units){ /*Get maximum:*/ - double maxvz=this->inputs->Max(VzEnum); + IssmDouble maxvz=this->inputs->Max(VzEnum); /*process units if requested: */ if(process_units) maxvz=UnitConversion(maxvz,IuToExtEnum,VzEnum); @@ -2796,10 +2796,10 @@ } /*}}}*/ /*FUNCTION Tria::MinVel{{{*/ -void Tria::MinVel(double* pminvel, bool process_units){ +void Tria::MinVel(IssmDouble* pminvel, bool process_units){ /*Get minimum:*/ - double minvel=this->inputs->Min(VelEnum); + IssmDouble minvel=this->inputs->Min(VelEnum); /*process units if requested: */ if(process_units) minvel=UnitConversion(minvel,IuToExtEnum,VelEnum); @@ -2809,10 +2809,10 @@ } /*}}}*/ /*FUNCTION Tria::MinVx{{{*/ -void Tria::MinVx(double* pminvx, bool process_units){ +void Tria::MinVx(IssmDouble* pminvx, bool process_units){ /*Get minimum:*/ - double minvx=this->inputs->Min(VxEnum); + IssmDouble minvx=this->inputs->Min(VxEnum); /*process units if requested: */ if(process_units) minvx=UnitConversion(minvx,IuToExtEnum,VxEnum); @@ -2822,10 +2822,10 @@ } /*}}}*/ /*FUNCTION Tria::MinVy{{{*/ -void Tria::MinVy(double* pminvy, bool process_units){ +void Tria::MinVy(IssmDouble* pminvy, bool process_units){ /*Get minimum:*/ - double minvy=this->inputs->Min(VyEnum); + IssmDouble minvy=this->inputs->Min(VyEnum); /*process units if requested: */ if(process_units) minvy=UnitConversion(minvy,IuToExtEnum,VyEnum); @@ -2835,10 +2835,10 @@ } /*}}}*/ /*FUNCTION Tria::MinVz{{{*/ -void Tria::MinVz(double* pminvz, bool process_units){ +void Tria::MinVz(IssmDouble* pminvz, bool process_units){ /*Get minimum:*/ - double minvz=this->inputs->Min(VzEnum); + IssmDouble minvz=this->inputs->Min(VzEnum); /*process units if requested: */ if(process_units) minvz=UnitConversion(minvz,IuToExtEnum,VzEnum); @@ -2848,7 +2848,7 @@ } /*}}}*/ /*FUNCTION Tria::ElementResponse{{{*/ -void Tria::ElementResponse(double* presponse,int response_enum,bool process_units){ +void Tria::ElementResponse(IssmDouble* presponse,int response_enum,bool process_units){ switch(response_enum){ case MaterialsRheologyBbarEnum: @@ -2857,7 +2857,7 @@ case VelEnum: /*Get input:*/ - double vel; + IssmDouble vel; Input* vel_input; vel_input=this->inputs->GetInput(VelEnum); _assert_(vel_input); @@ -2899,14 +2899,14 @@ /*Intermediaries*/ int i,j,ig; - double xyz_list[NUMVERTICES][3]; - double viscosity,newviscosity,oldviscosity; - double viscosity_overshoot,thickness,Jdet; - double epsilon[3],oldepsilon[3]; /* epsilon=[exx,eyy,exy]; */ - double B[3][numdof]; - double Bprime[3][numdof]; - double D[3][3] = {0.0}; - double D_scalar; + IssmDouble xyz_list[NUMVERTICES][3]; + IssmDouble viscosity,newviscosity,oldviscosity; + IssmDouble viscosity_overshoot,thickness,Jdet; + IssmDouble epsilon[3],oldepsilon[3]; /* epsilon=[exx,eyy,exy]; */ + IssmDouble B[3][numdof]; + IssmDouble Bprime[3][numdof]; + IssmDouble D[3][3] = {0.0}; + IssmDouble D_scalar; GaussTria *gauss = NULL; /*Initialize Element matrix*/ @@ -2964,15 +2964,15 @@ /*Intermediaries*/ int i,j,ig; int analysis_type; - double MAXSLOPE = .06; // 6 % - double MOUNTAINKEXPONENT = 10; - double slope_magnitude,alpha2; - double Jdet; - double L[2][numdof]; - double DL[2][2] = {{ 0,0 },{0,0}}; - double DL_scalar; - double slope[2] = {0.0,0.0}; - double xyz_list[NUMVERTICES][3]; + IssmDouble MAXSLOPE = .06; // 6 % + IssmDouble MOUNTAINKEXPONENT = 10; + IssmDouble slope_magnitude,alpha2; + IssmDouble Jdet; + IssmDouble L[2][numdof]; + IssmDouble DL[2][2] = {{ 0,0 },{0,0}}; + IssmDouble DL_scalar; + IssmDouble slope[2] = {0.0,0.0}; + IssmDouble xyz_list[NUMVERTICES][3]; Friction *friction = NULL; GaussTria *gauss = NULL; @@ -3001,7 +3001,7 @@ //velocity should be = 0. To achieve this result, we set alpha2_list to a very high value: */ surface_input->GetInputDerivativeValue(&slope[0],&xyz_list[0][0],gauss); slope_magnitude=sqrt(pow(slope[0],2)+pow(slope[1],2)); - if(slope_magnitude>MAXSLOPE) alpha2=pow((double)10,MOUNTAINKEXPONENT); + if(slope_magnitude>MAXSLOPE) alpha2=pow((IssmDouble)10,MOUNTAINKEXPONENT); else friction->GetAlpha2(&alpha2, gauss,VxEnum,VyEnum,VzEnum); GetL(&L[0][0], &xyz_list[0][0], gauss,NDOF2); @@ -3037,8 +3037,8 @@ /*Create Element matrix*/ for(i=0;iGetConnectivity(); - Ke->values[(2*i)*numdof +(2*i) ]=1/(double)connectivity; - Ke->values[(2*i+1)*numdof+(2*i+1)]=1/(double)connectivity; + Ke->values[(2*i)*numdof +(2*i) ]=1/(IssmDouble)connectivity; + Ke->values[(2*i+1)*numdof+(2*i+1)]=1/(IssmDouble)connectivity; } /*Clean up and return*/ @@ -3053,12 +3053,12 @@ /*Intermediaries */ int i,j,ig; - double driving_stress_baseline,thickness; - double Jdet; - double xyz_list[NUMVERTICES][3]; - double slope[2]; - double basis[3]; - double pe_g_gaussian[numdof]; + IssmDouble driving_stress_baseline,thickness; + IssmDouble Jdet; + IssmDouble xyz_list[NUMVERTICES][3]; + IssmDouble slope[2]; + IssmDouble basis[3]; + IssmDouble pe_g_gaussian[numdof]; GaussTria* gauss=NULL; /*Initialize Element vector*/ @@ -3104,10 +3104,10 @@ /*Intermediaries */ int i,connectivity; - double constant_part,ub,vb; - double rho_ice,gravity,n,B; - double slope2,thickness; - double slope[2]; + IssmDouble constant_part,ub,vb; + IssmDouble rho_ice,gravity,n,B; + IssmDouble slope2,thickness; + IssmDouble slope[2]; GaussTria* gauss=NULL; /*Initialize Element vector*/ @@ -3137,11 +3137,11 @@ constant_part=-2*pow(rho_ice*gravity,n)*pow(slope2,((n-1)/2)); - ub=-1.58*pow((double)10.0,(double)-10.0)*rho_ice*gravity*thickness*slope[0]; - vb=-1.58*pow((double)10.0,(double)-10.0)*rho_ice*gravity*thickness*slope[1]; + ub=-1.58*pow((IssmDouble)10.0,(IssmDouble)-10.0)*rho_ice*gravity*thickness*slope[0]; + vb=-1.58*pow((IssmDouble)10.0,(IssmDouble)-10.0)*rho_ice*gravity*thickness*slope[1]; - pe->values[2*i] =(ub-2.0*pow(rho_ice*gravity,n)*pow(slope2,((n-1)/2.0))*pow(thickness,n)/(pow(B,n)*(n+1))*slope[0])/(double)connectivity; - pe->values[2*i+1]=(vb-2.0*pow(rho_ice*gravity,n)*pow(slope2,((n-1)/2.0))*pow(thickness,n)/(pow(B,n)*(n+1))*slope[1])/(double)connectivity; + pe->values[2*i] =(ub-2.0*pow(rho_ice*gravity,n)*pow(slope2,((n-1)/2.0))*pow(thickness,n)/(pow(B,n)*(n+1))*slope[0])/(IssmDouble)connectivity; + pe->values[2*i+1]=(vb-2.0*pow(rho_ice*gravity,n)*pow(slope2,((n-1)/2.0))*pow(thickness,n)/(pow(B,n)*(n+1))*slope[1])/(IssmDouble)connectivity; } /*Clean up and return*/ @@ -3157,15 +3157,15 @@ /*Intermediaries */ int i,j,ig; - double xyz_list[NUMVERTICES][3]; - double Jdet,thickness; - double eps1dotdphii,eps1dotdphij; - double eps2dotdphii,eps2dotdphij; - double mu_prime; - double epsilon[3];/* epsilon=[exx,eyy,exy];*/ - double eps1[2],eps2[2]; - double phi[NUMVERTICES]; - double dphi[2][NUMVERTICES]; + IssmDouble xyz_list[NUMVERTICES][3]; + IssmDouble Jdet,thickness; + IssmDouble eps1dotdphii,eps1dotdphij; + IssmDouble eps2dotdphii,eps2dotdphij; + IssmDouble mu_prime; + IssmDouble epsilon[3];/* epsilon=[exx,eyy,exy];*/ + IssmDouble eps1[2],eps2[2]; + IssmDouble phi[NUMVERTICES]; + IssmDouble dphi[2][NUMVERTICES]; GaussTria *gauss=NULL; /*Initialize Jacobian with regular MacAyeal (first part of the Gateau derivative)*/ @@ -3222,8 +3222,8 @@ int i; int* doflist=NULL; - double vx,vy; - double values[numdof]; + IssmDouble vx,vy; + IssmDouble values[numdof]; GaussTria* gauss=NULL; /*Get dof list: */ @@ -3260,8 +3260,8 @@ const int numdof=NDOF2*NUMVERTICES; int i; - double vx,vy; - double values[numdof]; + IssmDouble vx,vy; + IssmDouble values[numdof]; int *doflist = NULL; GaussTria *gauss = NULL; @@ -3294,20 +3294,20 @@ } /*}}}*/ /*FUNCTION Tria::InputUpdateFromSolutionDiagnosticHoriz {{{*/ -void Tria::InputUpdateFromSolutionDiagnosticHoriz(double* solution){ +void Tria::InputUpdateFromSolutionDiagnosticHoriz(IssmDouble* solution){ const int numdof=NDOF2*NUMVERTICES; int i; int* doflist=NULL; - double rho_ice,g; - double values[numdof]; - double vx[NUMVERTICES]; - double vy[NUMVERTICES]; - double vz[NUMVERTICES]; - double vel[NUMVERTICES]; - double pressure[NUMVERTICES]; - double thickness[NUMVERTICES]; + IssmDouble rho_ice,g; + IssmDouble values[numdof]; + IssmDouble vx[NUMVERTICES]; + IssmDouble vy[NUMVERTICES]; + IssmDouble vz[NUMVERTICES]; + IssmDouble vel[NUMVERTICES]; + IssmDouble pressure[NUMVERTICES]; + IssmDouble thickness[NUMVERTICES]; /*Get dof list: */ GetDofList(&doflist,NoneApproximationEnum,GsetEnum); @@ -3357,20 +3357,20 @@ } /*}}}*/ /*FUNCTION Tria::InputUpdateFromSolutionDiagnosticHutter {{{*/ -void Tria::InputUpdateFromSolutionDiagnosticHutter(double* solution){ +void Tria::InputUpdateFromSolutionDiagnosticHutter(IssmDouble* solution){ const int numdof=NDOF2*NUMVERTICES; int i; int* doflist=NULL; - double rho_ice,g; - double values[numdof]; - double vx[NUMVERTICES]; - double vy[NUMVERTICES]; - double vz[NUMVERTICES]; - double vel[NUMVERTICES]; - double pressure[NUMVERTICES]; - double thickness[NUMVERTICES]; + IssmDouble rho_ice,g; + IssmDouble values[numdof]; + IssmDouble vx[NUMVERTICES]; + IssmDouble vy[NUMVERTICES]; + IssmDouble vz[NUMVERTICES]; + IssmDouble vel[NUMVERTICES]; + IssmDouble pressure[NUMVERTICES]; + IssmDouble thickness[NUMVERTICES]; /*Get dof list: */ GetDofList(&doflist,NoneApproximationEnum,GsetEnum); @@ -3419,7 +3419,7 @@ #ifdef _HAVE_CONTROL_ /*FUNCTION Tria::InputControlUpdate{{{*/ -void Tria::InputControlUpdate(double scalar,bool save_parameter){ +void Tria::InputControlUpdate(IssmDouble scalar,bool save_parameter){ /*Intermediary*/ int num_controls; @@ -3473,7 +3473,7 @@ }/*}}}*/ /*FUNCTION Tria::ControlInputScaleGradient{{{*/ -void Tria::ControlInputScaleGradient(int enum_type,double scale){ +void Tria::ControlInputScaleGradient(int enum_type,IssmDouble scale){ Input* input=NULL; @@ -3489,10 +3489,10 @@ ((ControlInput*)input)->ScaleGradient(scale); }/*}}}*/ /*FUNCTION Tria::ControlInputSetGradient{{{*/ -void Tria::ControlInputSetGradient(double* gradient,int enum_type,int control_index){ +void Tria::ControlInputSetGradient(IssmDouble* gradient,int enum_type,int control_index){ int doflist1[NUMVERTICES]; - double grad_list[NUMVERTICES]; + IssmDouble grad_list[NUMVERTICES]; Input* grad_input=NULL; Input* input=NULL; @@ -3576,11 +3576,11 @@ int i,ig; int doflist1[NUMVERTICES]; - double Jdet,weight; - double xyz_list[NUMVERTICES][3]; - double dbasis[NDOF2][NUMVERTICES]; - double dk[NDOF2]; - double grade_g[NUMVERTICES]={0.0}; + IssmDouble Jdet,weight; + IssmDouble xyz_list[NUMVERTICES][3]; + IssmDouble dbasis[NDOF2][NUMVERTICES]; + IssmDouble dk[NDOF2]; + IssmDouble grade_g[NUMVERTICES]={0.0}; GaussTria *gauss=NULL; /*Retrieve all inputs we will be needing: */ @@ -3617,12 +3617,12 @@ /*Intermediaries*/ int i,ig; int doflist[NUMVERTICES]; - double vx,vy,lambda,mu,thickness,Jdet; - double viscosity_complement; - double dvx[NDOF2],dvy[NDOF2],dadjx[NDOF2],dadjy[NDOF2],dB[NDOF2]; - double xyz_list[NUMVERTICES][3]; - double basis[3],epsilon[3]; - double grad[NUMVERTICES]={0.0}; + IssmDouble vx,vy,lambda,mu,thickness,Jdet; + IssmDouble viscosity_complement; + IssmDouble dvx[NDOF2],dvy[NDOF2],dadjx[NDOF2],dadjy[NDOF2],dB[NDOF2]; + IssmDouble xyz_list[NUMVERTICES][3]; + IssmDouble basis[3],epsilon[3]; + IssmDouble grad[NUMVERTICES]={0.0}; GaussTria *gauss = NULL; /* Get node coordinates and dof list: */ @@ -3675,14 +3675,14 @@ int analysis_type; int doflist1[NUMVERTICES]; int connectivity[NUMVERTICES]; - double vx,vy,lambda,mu,alpha_complement,Jdet; - double bed,thickness,Neff,drag; - double xyz_list[NUMVERTICES][3]; - double dk[NDOF2]; - double grade_g[NUMVERTICES]={0.0}; - double grade_g_gaussian[NUMVERTICES]; - double basis[3]; - double epsilon[3]; /* epsilon=[exx,eyy,exy];*/ + IssmDouble vx,vy,lambda,mu,alpha_complement,Jdet; + IssmDouble bed,thickness,Neff,drag; + IssmDouble xyz_list[NUMVERTICES][3]; + IssmDouble dk[NDOF2]; + IssmDouble grade_g[NUMVERTICES]={0.0}; + IssmDouble grade_g_gaussian[NUMVERTICES]; + IssmDouble basis[3]; + IssmDouble epsilon[3]; /* epsilon=[exx,eyy,exy];*/ Friction* friction=NULL; GaussTria *gauss=NULL; @@ -3745,7 +3745,7 @@ // adjointy_input->GetInputValue(&mu, gauss); // vx_input->GetInputValue(&vx,gauss); // vy_input->GetInputValue(&vy,gauss); - // grade_g[iv] = -2*1.e+7*drag*alpha_complement*(lambda*vx+mu*vy)/((double)connectivity[iv]); + // grade_g[iv] = -2*1.e+7*drag*alpha_complement*(lambda*vx+mu*vy)/((IssmDouble)connectivity[iv]); //} /*End Analytical gradient*/ @@ -3761,11 +3761,11 @@ int i,ig; int doflist1[NUMVERTICES]; - double Jdet,weight; - double xyz_list[NUMVERTICES][3]; - double dbasis[NDOF2][NUMVERTICES]; - double dk[NDOF2]; - double grade_g[NUMVERTICES]={0.0}; + IssmDouble Jdet,weight; + IssmDouble xyz_list[NUMVERTICES][3]; + IssmDouble dbasis[NDOF2][NUMVERTICES]; + IssmDouble dk[NDOF2]; + IssmDouble grade_g[NUMVERTICES]={0.0}; GaussTria *gauss=NULL; /*Retrieve all inputs we will be needing: */ @@ -3805,8 +3805,8 @@ /*Intermediaries*/ int doflist1[NUMVERTICES]; - double lambda[NUMVERTICES]; - double gradient_g[NUMVERTICES]; + IssmDouble lambda[NUMVERTICES]; + IssmDouble gradient_g[NUMVERTICES]; /*Compute Gradient*/ GradientIndexing(&doflist1[0],control_index); @@ -3822,11 +3822,11 @@ /*Intermediaries*/ int i,ig; int doflist1[NUMVERTICES]; - double thickness,Jdet; - double basis[3]; - double Dlambda[2],dp[2]; - double xyz_list[NUMVERTICES][3]; - double grade_g[NUMVERTICES] = {0.0}; + IssmDouble thickness,Jdet; + IssmDouble basis[3]; + IssmDouble Dlambda[2],dp[2]; + IssmDouble xyz_list[NUMVERTICES][3]; + IssmDouble grade_g[NUMVERTICES] = {0.0}; GaussTria *gauss = NULL; /* Get node coordinates and dof list: */ @@ -3865,11 +3865,11 @@ /*Intermediaries*/ int i,ig; int doflist1[NUMVERTICES]; - double thickness,Jdet; - double basis[3]; - double Dlambda[2],dp[2]; - double xyz_list[NUMVERTICES][3]; - double grade_g[NUMVERTICES] = {0.0}; + IssmDouble thickness,Jdet; + IssmDouble basis[3]; + IssmDouble Dlambda[2],dp[2]; + IssmDouble xyz_list[NUMVERTICES][3]; + IssmDouble grade_g[NUMVERTICES] = {0.0}; GaussTria *gauss = NULL; /* Get node coordinates and dof list: */ @@ -3916,15 +3916,15 @@ } /*}}}*/ /*FUNCTION Tria::RheologyBbarAbsGradient{{{*/ -double Tria::RheologyBbarAbsGradient(bool process_units,int weight_index){ +IssmDouble Tria::RheologyBbarAbsGradient(bool process_units,int weight_index){ /* Intermediaries */ int ig; - double Jelem = 0; - double weight; - double Jdet; - double xyz_list[NUMVERTICES][3]; - double dp[NDOF2]; + IssmDouble Jelem = 0; + IssmDouble weight; + IssmDouble Jdet; + IssmDouble xyz_list[NUMVERTICES][3]; + IssmDouble dp[NDOF2]; GaussTria *gauss = NULL; /*retrieve parameters and inputs*/ @@ -3960,15 +3960,15 @@ } /*}}}*/ /*FUNCTION Tria::SurfaceAverageVelMisfit {{{*/ -double Tria::SurfaceAverageVelMisfit(bool process_units,int weight_index){ +IssmDouble Tria::SurfaceAverageVelMisfit(bool process_units,int weight_index){ const int numdof=2*NUMVERTICES; int i,ig; - double Jelem=0,S,Jdet; - double misfit; - double vx,vy,vxobs,vyobs,weight; - double xyz_list[NUMVERTICES][3]; + IssmDouble Jelem=0,S,Jdet; + IssmDouble misfit; + IssmDouble vx,vy,vxobs,vyobs,weight; + IssmDouble xyz_list[NUMVERTICES][3]; GaussTria *gauss=NULL; /*If on water, return 0: */ @@ -4021,18 +4021,18 @@ } /*}}}*/ /*FUNCTION Tria::SurfaceLogVelMisfit {{{*/ -double Tria::SurfaceLogVelMisfit(bool process_units,int weight_index){ +IssmDouble Tria::SurfaceLogVelMisfit(bool process_units,int weight_index){ const int numdof=NDOF2*NUMVERTICES; int i,ig; - double Jelem=0; - double misfit,Jdet; - double epsvel=2.220446049250313e-16; - double meanvel=3.170979198376458e-05; /*1000 m/yr*/ - double velocity_mag,obs_velocity_mag; - double xyz_list[NUMVERTICES][3]; - double vx,vy,vxobs,vyobs,weight; + IssmDouble Jelem=0; + IssmDouble misfit,Jdet; + IssmDouble epsvel=2.220446049250313e-16; + IssmDouble meanvel=3.170979198376458e-05; /*1000 m/yr*/ + IssmDouble velocity_mag,obs_velocity_mag; + IssmDouble xyz_list[NUMVERTICES][3]; + IssmDouble vx,vy,vxobs,vyobs,weight; GaussTria *gauss=NULL; /*If on water, return 0: */ @@ -4086,18 +4086,18 @@ } /*}}}*/ /*FUNCTION Tria::SurfaceLogVxVyMisfit {{{*/ -double Tria::SurfaceLogVxVyMisfit(bool process_units,int weight_index){ +IssmDouble Tria::SurfaceLogVxVyMisfit(bool process_units,int weight_index){ const int numdof=NDOF2*NUMVERTICES; int i,ig; int fit=-1; - double Jelem=0, S=0; - double epsvel=2.220446049250313e-16; - double meanvel=3.170979198376458e-05; /*1000 m/yr*/ - double misfit, Jdet; - double vx,vy,vxobs,vyobs,weight; - double xyz_list[NUMVERTICES][3]; + IssmDouble Jelem=0, S=0; + IssmDouble epsvel=2.220446049250313e-16; + IssmDouble meanvel=3.170979198376458e-05; /*1000 m/yr*/ + IssmDouble misfit, Jdet; + IssmDouble vx,vy,vxobs,vyobs,weight; + IssmDouble xyz_list[NUMVERTICES][3]; GaussTria *gauss=NULL; /*If on water, return 0: */ @@ -4152,15 +4152,15 @@ } /*}}}*/ /*FUNCTION Tria::SurfaceAbsVelMisfit {{{*/ -double Tria::SurfaceAbsVelMisfit(bool process_units,int weight_index){ +IssmDouble Tria::SurfaceAbsVelMisfit(bool process_units,int weight_index){ const int numdof=NDOF2*NUMVERTICES; int i,ig; - double Jelem=0; - double misfit,Jdet; - double vx,vy,vxobs,vyobs,weight; - double xyz_list[NUMVERTICES][3]; + IssmDouble Jelem=0; + IssmDouble misfit,Jdet; + IssmDouble vx,vy,vxobs,vyobs,weight; + IssmDouble xyz_list[NUMVERTICES][3]; GaussTria *gauss=NULL; /*If on water, return 0: */ @@ -4213,17 +4213,17 @@ } /*}}}*/ /*FUNCTION Tria::SurfaceRelVelMisfit {{{*/ -double Tria::SurfaceRelVelMisfit(bool process_units,int weight_index){ +IssmDouble Tria::SurfaceRelVelMisfit(bool process_units,int weight_index){ const int numdof=2*NUMVERTICES; int i,ig; - double Jelem=0; - double scalex=1,scaley=1; - double misfit,Jdet; - double epsvel=2.220446049250313e-16; - double meanvel=3.170979198376458e-05; /*1000 m/yr*/ - double vx,vy,vxobs,vyobs,weight; - double xyz_list[NUMVERTICES][3]; + IssmDouble Jelem=0; + IssmDouble scalex=1,scaley=1; + IssmDouble misfit,Jdet; + IssmDouble epsvel=2.220446049250313e-16; + IssmDouble meanvel=3.170979198376458e-05; /*1000 m/yr*/ + IssmDouble vx,vy,vxobs,vyobs,weight; + IssmDouble xyz_list[NUMVERTICES][3]; GaussTria *gauss=NULL; /*If on water, return 0: */ @@ -4277,15 +4277,15 @@ } /*}}}*/ /*FUNCTION Tria::ThicknessAbsGradient{{{*/ -double Tria::ThicknessAbsGradient(bool process_units,int weight_index){ +IssmDouble Tria::ThicknessAbsGradient(bool process_units,int weight_index){ /* Intermediaries */ int ig; - double Jelem = 0; - double weight; - double Jdet; - double xyz_list[NUMVERTICES][3]; - double dp[NDOF2]; + IssmDouble Jelem = 0; + IssmDouble weight; + IssmDouble Jdet; + IssmDouble xyz_list[NUMVERTICES][3]; + IssmDouble dp[NDOF2]; GaussTria *gauss = NULL; /*retrieve parameters and inputs*/ @@ -4321,16 +4321,16 @@ } /*}}}*/ /*FUNCTION Tria::ThicknessAbsMisfit {{{*/ -double Tria::ThicknessAbsMisfit(bool process_units,int weight_index){ +IssmDouble Tria::ThicknessAbsMisfit(bool process_units,int weight_index){ /*Intermediaries*/ int i,ig; - double thickness,thicknessobs,weight; - double Jdet; - double Jelem = 0; - double xyz_list[NUMVERTICES][3]; + IssmDouble thickness,thicknessobs,weight; + IssmDouble Jdet; + IssmDouble Jelem = 0; + IssmDouble xyz_list[NUMVERTICES][3]; GaussTria *gauss = NULL; - double dH[2]; + IssmDouble dH[2]; /*If on water, return 0: */ if(IsOnWater())return 0; @@ -4373,14 +4373,14 @@ /*Intermediaries */ int i,ig,resp; - double Jdet; - double thickness,thicknessobs,weight; + IssmDouble Jdet; + IssmDouble thickness,thicknessobs,weight; int *responses = NULL; int num_responses; - double xyz_list[NUMVERTICES][3]; - double basis[3]; - double dbasis[NDOF2][NUMVERTICES]; - double dH[2]; + IssmDouble xyz_list[NUMVERTICES][3]; + IssmDouble basis[3]; + IssmDouble dbasis[NDOF2][NUMVERTICES]; + IssmDouble dH[2]; GaussTria* gauss=NULL; /*Initialize Element vector*/ @@ -4441,15 +4441,15 @@ int i,resp,ig; int *responses=NULL; int num_responses; - double Jdet; - double obs_velocity_mag,velocity_mag; - double dux,duy; - double epsvel=2.220446049250313e-16; - double meanvel=3.170979198376458e-05; /*1000 m/yr*/ - double scalex=0,scaley=0,scale=0,S=0; - double vx,vy,vxobs,vyobs,weight; - double xyz_list[NUMVERTICES][3]; - double basis[3]; + IssmDouble Jdet; + IssmDouble obs_velocity_mag,velocity_mag; + IssmDouble dux,duy; + IssmDouble epsvel=2.220446049250313e-16; + IssmDouble meanvel=3.170979198376458e-05; /*1000 m/yr*/ + IssmDouble scalex=0,scaley=0,scale=0,S=0; + IssmDouble vx,vy,vxobs,vyobs,weight; + IssmDouble xyz_list[NUMVERTICES][3]; + IssmDouble basis[3]; GaussTria* gauss=NULL; /*Initialize Element vector*/ @@ -4617,15 +4617,15 @@ int i,resp,ig; int *responses=NULL; int num_responses; - double Jdet; - double obs_velocity_mag,velocity_mag; - double dux,duy; - double epsvel=2.220446049250313e-16; - double meanvel=3.170979198376458e-05; /*1000 m/yr*/ - double scalex=0,scaley=0,scale=0,S=0; - double vx,vy,vxobs,vyobs,weight; - double xyz_list[NUMVERTICES][3]; - double basis[3]; + IssmDouble Jdet; + IssmDouble obs_velocity_mag,velocity_mag; + IssmDouble dux,duy; + IssmDouble epsvel=2.220446049250313e-16; + IssmDouble meanvel=3.170979198376458e-05; /*1000 m/yr*/ + IssmDouble scalex=0,scaley=0,scale=0,S=0; + IssmDouble vx,vy,vxobs,vyobs,weight; + IssmDouble xyz_list[NUMVERTICES][3]; + IssmDouble basis[3]; GaussTria* gauss=NULL; /*Initialize Element vector*/ @@ -4788,15 +4788,15 @@ } /*}}}*/ /*FUNCTION Tria::DragCoefficientAbsGradient{{{*/ -double Tria::DragCoefficientAbsGradient(bool process_units,int weight_index){ +IssmDouble Tria::DragCoefficientAbsGradient(bool process_units,int weight_index){ /* Intermediaries */ int ig; - double Jelem = 0; - double weight; - double Jdet; - double xyz_list[NUMVERTICES][3]; - double dp[NDOF2]; + IssmDouble Jelem = 0; + IssmDouble weight; + IssmDouble Jdet; + IssmDouble xyz_list[NUMVERTICES][3]; + IssmDouble dp[NDOF2]; GaussTria *gauss = NULL; /*retrieve parameters and inputs*/ @@ -4862,15 +4862,15 @@ /*Intermediaries */ int i,j,ig; bool incomplete_adjoint; - double xyz_list[NUMVERTICES][3]; - double Jdet,thickness; - double eps1dotdphii,eps1dotdphij; - double eps2dotdphii,eps2dotdphij; - double mu_prime; - double epsilon[3];/* epsilon=[exx,eyy,exy];*/ - double eps1[2],eps2[2]; - double phi[NUMVERTICES]; - double dphi[2][NUMVERTICES]; + IssmDouble xyz_list[NUMVERTICES][3]; + IssmDouble Jdet,thickness; + IssmDouble eps1dotdphii,eps1dotdphij; + IssmDouble eps2dotdphii,eps2dotdphij; + IssmDouble mu_prime; + IssmDouble epsilon[3];/* epsilon=[exx,eyy,exy];*/ + IssmDouble eps1[2],eps2[2]; + IssmDouble phi[NUMVERTICES]; + IssmDouble dphi[2][NUMVERTICES]; GaussTria *gauss=NULL; /*Initialize Jacobian with regular MacAyeal (first part of the Gateau derivative)*/ @@ -4924,15 +4924,15 @@ } /*}}}*/ /*FUNCTION Tria::InputUpdateFromSolutionAdjointHoriz {{{*/ -void Tria::InputUpdateFromSolutionAdjointHoriz(double* solution){ +void Tria::InputUpdateFromSolutionAdjointHoriz(IssmDouble* solution){ const int numdof=NDOF2*NUMVERTICES; int i; int* doflist=NULL; - double values[numdof]; - double lambdax[NUMVERTICES]; - double lambday[NUMVERTICES]; + IssmDouble values[numdof]; + IssmDouble lambdax[NUMVERTICES]; + IssmDouble lambday[NUMVERTICES]; /*Get dof list: */ GetDofList(&doflist,NoneApproximationEnum,GsetEnum); @@ -4959,14 +4959,14 @@ } /*}}}*/ /*FUNCTION Tria::InputUpdateFromSolutionAdjointBalancethickness {{{*/ -void Tria::InputUpdateFromSolutionAdjointBalancethickness(double* solution){ +void Tria::InputUpdateFromSolutionAdjointBalancethickness(IssmDouble* solution){ const int numdof=NDOF1*NUMVERTICES; int i; int* doflist=NULL; - double values[numdof]; - double lambda[NUMVERTICES]; + IssmDouble values[numdof]; + IssmDouble lambda[NUMVERTICES]; /*Get dof list: */ GetDofList(&doflist,NoneApproximationEnum,GsetEnum); @@ -5016,9 +5016,9 @@ } /*}}}*/ /*FUNCTION Tria::SetControlInputsFromVector{{{*/ -void Tria::SetControlInputsFromVector(double* vector,int control_enum,int control_index){ +void Tria::SetControlInputsFromVector(IssmDouble* vector,int control_enum,int control_index){ - double values[NUMVERTICES]; + IssmDouble values[NUMVERTICES]; int doflist1[NUMVERTICES]; Input *input = NULL; Input *new_input = NULL; @@ -5056,14 +5056,14 @@ void Tria::CreateHydrologyWaterVelocityInput(void){ /*material parameters: */ - double mu_water; - double VelocityFactor; // This factor represents the number 12 in laminar flow velocity which can vary by differnt hydrology.CR - double n_man,CR; - double w; - double rho_ice, rho_water, g; - double dsdx,dsdy,dbdx,dbdy; - double vx[NUMVERTICES]; - double vy[NUMVERTICES]; + IssmDouble mu_water; + IssmDouble VelocityFactor; // This factor represents the number 12 in laminar flow velocity which can vary by differnt hydrology.CR + IssmDouble n_man,CR; + IssmDouble w; + IssmDouble rho_ice, rho_water, g; + IssmDouble dsdx,dsdy,dbdx,dbdy; + IssmDouble vx[NUMVERTICES]; + IssmDouble vy[NUMVERTICES]; GaussTria *gauss = NULL; /*Retrieve all inputs and parameters*/ @@ -5114,20 +5114,20 @@ const int numdof=NDOF1*NUMVERTICES; /*Intermediaries */ - double diffusivity; + IssmDouble diffusivity; int i,j,ig; - double Jdettria,DL_scalar,dt,h; - double vx,vy,vel,dvxdx,dvydy; - double dvx[2],dvy[2]; - double v_gauss[2]={0.0}; - double xyz_list[NUMVERTICES][3]; - double L[NUMVERTICES]; - double B[2][NUMVERTICES]; - double Bprime[2][NUMVERTICES]; - double K[2][2] ={0.0}; - double KDL[2][2] ={0.0}; - double DL[2][2] ={0.0}; - double DLprime[2][2] ={0.0}; + IssmDouble Jdettria,DL_scalar,dt,h; + IssmDouble vx,vy,vel,dvxdx,dvydy; + IssmDouble dvx[2],dvy[2]; + IssmDouble v_gauss[2]={0.0}; + IssmDouble xyz_list[NUMVERTICES][3]; + IssmDouble L[NUMVERTICES]; + IssmDouble B[2][NUMVERTICES]; + IssmDouble Bprime[2][NUMVERTICES]; + IssmDouble K[2][2] ={0.0}; + IssmDouble KDL[2][2] ={0.0}; + IssmDouble DL[2][2] ={0.0}; + IssmDouble DLprime[2][2] ={0.0}; GaussTria *gauss=NULL; /*Skip if water or ice shelf element*/ @@ -5220,11 +5220,11 @@ /*Intermediaries */ int i,j,ig; - double Jdettria,dt; - double basal_melting_g; - double old_watercolumn_g; - double xyz_list[NUMVERTICES][3]; - double basis[numdof]; + IssmDouble Jdettria,dt; + IssmDouble basal_melting_g; + IssmDouble old_watercolumn_g; + IssmDouble xyz_list[NUMVERTICES][3]; + IssmDouble basis[numdof]; GaussTria* gauss=NULL; /*Skip if water or ice shelf element*/ @@ -5268,8 +5268,8 @@ int i; int* doflist=NULL; - double watercolumn; - double values[numdof]; + IssmDouble watercolumn; + IssmDouble values[numdof]; GaussTria* gauss=NULL; /*Get dof list: */ @@ -5298,14 +5298,14 @@ } /*}}}*/ /*FUNCTION Tria::InputUpdateFromSolutionHydrology{{{*/ -void Tria::InputUpdateFromSolutionHydrology(double* solution){ +void Tria::InputUpdateFromSolutionHydrology(IssmDouble* solution){ /*Intermediaries*/ const int numdof = NDOF1*NUMVERTICES; int i; int* doflist=NULL; - double values[numdof]; + IssmDouble values[numdof]; /*Get dof list: */ GetDofList(&doflist,NoneApproximationEnum,GsetEnum); @@ -5314,7 +5314,7 @@ for(i=0;iIsFloating()){ /*hydrostatic equilibrium: */ - double rho_ice,rho_water,di; + IssmDouble rho_ice,rho_water,di; rho_ice=this->matpar->GetRhoIce(); rho_water=this->matpar->GetRhoWater(); @@ -5440,15 +5440,15 @@ _error_(" not supported yet!"); } /*}}}*/ -/*FUNCTION Tria::InputUpdateFromMatrixDakota(double* matrix, int nrows, int ncols, int name, int type);{{{*/ -void Tria::InputUpdateFromMatrixDakota(double* matrix, int nrows, int ncols, int name, int type){ +/*FUNCTION Tria::InputUpdateFromMatrixDakota(IssmDouble* matrix, int nrows, int ncols, int name, int type);{{{*/ +void Tria::InputUpdateFromMatrixDakota(IssmDouble* matrix, int nrows, int ncols, int name, int type){ int i,j,t; TransientInput* transientinput=NULL; - double values[3]; - double time; + IssmDouble values[3]; + IssmDouble time; int row; - double yts; + IssmDouble yts; /*Check that name is an element input*/ if (!IsInput(name)) return; @@ -5465,11 +5465,11 @@ /*create input values: */ for(i=0;i<3;i++){ row=this->nodes[i]->GetSidList(); - values[i]=(double)matrix[ncols*row+t]; + values[i]=(IssmDouble)matrix[ncols*row+t]; } /*time? :*/ - time=(double)matrix[(nrows-1)*ncols+t]*yts; + time=(IssmDouble)matrix[(nrows-1)*ncols+t]*yts; if(t==0) transientinput=new TransientInput(name); transientinput->AddTimeInput(new TriaP1Input(name,values),time); @@ -5510,17 +5510,17 @@ /*Intermediaries */ int stabilization; int i,j,ig,dim; - double Jdettria,vx,vy,dvxdx,dvydy,vel,h; - double dvx[2],dvy[2]; - double xyz_list[NUMVERTICES][3]; - double L[NUMVERTICES]; - double B[2][NUMVERTICES]; - double Bprime[2][NUMVERTICES]; - double K[2][2] = {0.0}; - double KDL[2][2] = {0.0}; - double DL[2][2] = {0.0}; - double DLprime[2][2] = {0.0}; - double DL_scalar; + IssmDouble Jdettria,vx,vy,dvxdx,dvydy,vel,h; + IssmDouble dvx[2],dvy[2]; + IssmDouble xyz_list[NUMVERTICES][3]; + IssmDouble L[NUMVERTICES]; + IssmDouble B[2][NUMVERTICES]; + IssmDouble Bprime[2][NUMVERTICES]; + IssmDouble K[2][2] = {0.0}; + IssmDouble KDL[2][2] = {0.0}; + IssmDouble DL[2][2] = {0.0}; + IssmDouble DLprime[2][2] = {0.0}; + IssmDouble DL_scalar; GaussTria *gauss = NULL; /*Initialize Element matrix*/ @@ -5619,12 +5619,12 @@ /*Intermediaries*/ int i,j,ig,dim; - double vx,vy,Jdettria; - double xyz_list[NUMVERTICES][3]; - double B[2][NUMVERTICES]; - double Bprime[2][NUMVERTICES]; - double DL[2][2]={0.0}; - double DL_scalar; + IssmDouble vx,vy,Jdettria; + IssmDouble xyz_list[NUMVERTICES][3]; + IssmDouble B[2][NUMVERTICES]; + IssmDouble Bprime[2][NUMVERTICES]; + IssmDouble DL[2][2]={0.0}; + IssmDouble DL_scalar; GaussTria *gauss=NULL; /*Initialize Element matrix*/ @@ -5687,9 +5687,9 @@ /*Intermediaries */ int i,j,ig; - double xyz_list[NUMVERTICES][3]; - double dhdt_g,basal_melting_g,surface_mass_balance_g,Jdettria; - double L[NUMVERTICES]; + IssmDouble xyz_list[NUMVERTICES][3]; + IssmDouble dhdt_g,basal_melting_g,surface_mass_balance_g,Jdettria; + IssmDouble L[NUMVERTICES]; GaussTria* gauss=NULL; /*Initialize Element vector*/ @@ -5730,9 +5730,9 @@ /*Intermediaries */ int i,j,ig; - double xyz_list[NUMVERTICES][3]; - double basal_melting_g,surface_mass_balance_g,dhdt_g,Jdettria; - double L[NUMVERTICES]; + IssmDouble xyz_list[NUMVERTICES][3]; + IssmDouble basal_melting_g,surface_mass_balance_g,dhdt_g,Jdettria; + IssmDouble L[NUMVERTICES]; GaussTria* gauss=NULL; /*Initialize Element vector*/ Index: /u/astrid-r1b/morlighe/issmuci/trunk-jpl/../trunk-jpl/src/c/objects/Elements/TriaRef.h =================================================================== --- /u/astrid-r1b/morlighe/issmuci/trunk-jpl/../trunk-jpl/src/c/objects/Elements/TriaRef.h (revision 12470) +++ /u/astrid-r1b/morlighe/issmuci/trunk-jpl/../trunk-jpl/src/c/objects/Elements/TriaRef.h (revision 12471) @@ -24,26 +24,26 @@ void SetElementType(int type,int type_counter); /*Numerics*/ - void GetBMacAyeal(double* B, double* xyz_list, GaussTria* gauss); - void GetBMacAyealStokes(double* B , double* xyz_list, GaussTria* gauss); - void GetBprimeMacAyeal(double* Bprime, double* xyz_list, GaussTria* gauss); - void GetBprimeMacAyealStokes(double* Bprime, double* xyz_list, GaussTria* gauss); - void GetBprimePrognostic(double* Bprime_prog, double* xyz_list, GaussTria* gauss); - void GetBPrognostic(double* B_prog, double* xyz_list, GaussTria* gauss); - void GetL(double* L, double* xyz_list,GaussTria* gauss,int numdof); - void GetJacobian(double* J, double* xyz_list,GaussTria* gauss); - void GetSegmentJacobianDeterminant(double* Jdet, double* xyz_list,GaussTria* gauss); - void GetJacobianDeterminant2d(double* Jdet, double* xyz_list,GaussTria* gauss); - void GetJacobianDeterminant3d(double* Jdet, double* xyz_list,GaussTria* gauss); - void GetJacobianInvert(double* Jinv, double* xyz_list,GaussTria* gauss); - void GetNodalFunctions(double* l1l2l3,GaussTria* gauss); - void GetSegmentNodalFunctions(double* l1l2l3,GaussTria* gauss, int index1,int index2); - void GetSegmentBFlux(double* B,GaussTria* gauss, int index1,int index2); - void GetSegmentBprimeFlux(double* Bprime,GaussTria* gauss, int index1,int index2); - void GetNodalFunctionsDerivatives(double* l1l2l3,double* xyz_list, GaussTria* gauss); - void GetNodalFunctionsDerivativesReference(double* dl1dl3,GaussTria* gauss); - void GetInputValue(double* pp, double* plist, GaussTria* gauss); - void GetInputDerivativeValue(double* pp, double* plist,double* xyz_list, GaussTria* gauss); + void GetBMacAyeal(IssmDouble* B, IssmDouble* xyz_list, GaussTria* gauss); + void GetBMacAyealStokes(IssmDouble* B , IssmDouble* xyz_list, GaussTria* gauss); + void GetBprimeMacAyeal(IssmDouble* Bprime, IssmDouble* xyz_list, GaussTria* gauss); + void GetBprimeMacAyealStokes(IssmDouble* Bprime, IssmDouble* xyz_list, GaussTria* gauss); + void GetBprimePrognostic(IssmDouble* Bprime_prog, IssmDouble* xyz_list, GaussTria* gauss); + void GetBPrognostic(IssmDouble* B_prog, IssmDouble* xyz_list, GaussTria* gauss); + void GetL(IssmDouble* L, IssmDouble* xyz_list,GaussTria* gauss,int numdof); + void GetJacobian(IssmDouble* J, IssmDouble* xyz_list,GaussTria* gauss); + void GetSegmentJacobianDeterminant(IssmDouble* Jdet, IssmDouble* xyz_list,GaussTria* gauss); + void GetJacobianDeterminant2d(IssmDouble* Jdet, IssmDouble* xyz_list,GaussTria* gauss); + void GetJacobianDeterminant3d(IssmDouble* Jdet, IssmDouble* xyz_list,GaussTria* gauss); + void GetJacobianInvert(IssmDouble* Jinv, IssmDouble* xyz_list,GaussTria* gauss); + void GetNodalFunctions(IssmDouble* l1l2l3,GaussTria* gauss); + void GetSegmentNodalFunctions(IssmDouble* l1l2l3,GaussTria* gauss, int index1,int index2); + void GetSegmentBFlux(IssmDouble* B,GaussTria* gauss, int index1,int index2); + void GetSegmentBprimeFlux(IssmDouble* Bprime,GaussTria* gauss, int index1,int index2); + void GetNodalFunctionsDerivatives(IssmDouble* l1l2l3,IssmDouble* xyz_list, GaussTria* gauss); + void GetNodalFunctionsDerivativesReference(IssmDouble* dl1dl3,GaussTria* gauss); + void GetInputValue(IssmDouble* pp, IssmDouble* plist, GaussTria* gauss); + void GetInputDerivativeValue(IssmDouble* pp, IssmDouble* plist,IssmDouble* xyz_list, GaussTria* gauss); }; #endif Index: /u/astrid-r1b/morlighe/issmuci/trunk-jpl/../trunk-jpl/src/c/objects/Elements/PentaRef.cpp =================================================================== --- /u/astrid-r1b/morlighe/issmuci/trunk-jpl/../trunk-jpl/src/c/objects/Elements/PentaRef.cpp (revision 12470) +++ /u/astrid-r1b/morlighe/issmuci/trunk-jpl/../trunk-jpl/src/c/objects/Elements/PentaRef.cpp (revision 12471) @@ -57,7 +57,7 @@ /*Reference Element numerics*/ /*FUNCTION PentaRef::GetBMacAyealPattyn {{{*/ -void PentaRef::GetBMacAyealPattyn(double* B, double* xyz_list, GaussPenta* gauss){ +void PentaRef::GetBMacAyealPattyn(IssmDouble* B, IssmDouble* xyz_list, GaussPenta* gauss){ /*Compute B matrix. B=[B1 B2 B3 B4 B5 B6] where Bi is of size 5*NDOF2. * For node i, Bi can be expressed in the actual coordinate system * by: @@ -69,7 +69,7 @@ * We assume B has been allocated already, of size: 5x(NDOF2*NUMNODESP1) */ - double dbasis[3][NUMNODESP1]; + IssmDouble dbasis[3][NUMNODESP1]; /*Get dbasis in actual coordinate system: */ GetNodalFunctionsP1Derivatives(&dbasis[0][0],xyz_list, gauss); @@ -88,7 +88,7 @@ } /*}}}*/ /*FUNCTION PentaRef::GetBMacAyealStokes{{{*/ -void PentaRef::GetBMacAyealStokes(double* B, double* xyz_list, GaussPenta* gauss){ +void PentaRef::GetBMacAyealStokes(IssmDouble* B, IssmDouble* xyz_list, GaussPenta* gauss){ /*Compute B matrix. B=[B1 B2 B3 B4 B5 B6] where Bi is of size 5*NDOF2. * For node i, Bi can be expressed in the actual coordinate system * by: @@ -102,8 +102,8 @@ */ int i; - double dh1dh7[3][NUMNODESMINI]; - double l1l6[NUMNODESP1]; + IssmDouble dh1dh7[3][NUMNODESMINI]; + IssmDouble l1l6[NUMNODESP1]; /*Get dh1dh6 in actual coordinate system: */ GetNodalFunctionsMINIDerivatives(&dh1dh7[0][0],xyz_list, gauss); @@ -134,7 +134,7 @@ } /*}}}*/ /*FUNCTION PentaRef::GetBPattyn {{{*/ -void PentaRef::GetBPattyn(double* B, double* xyz_list, GaussPenta* gauss){ +void PentaRef::GetBPattyn(IssmDouble* B, IssmDouble* xyz_list, GaussPenta* gauss){ /*Compute B matrix. B=[B1 B2 B3 B4 B5 B6] where Bi is of size 5*NDOF2. * For node i, Bi can be expressed in the actual coordinate system * by: @@ -148,7 +148,7 @@ * We assume B has been allocated already, of size: 5x(NDOF2*NUMNODESP1) */ - double dbasis[3][NUMNODESP1]; + IssmDouble dbasis[3][NUMNODESP1]; /*Get dbasis in actual coordinate system: */ GetNodalFunctionsP1Derivatives(&dbasis[0][0],xyz_list, gauss); @@ -174,7 +174,7 @@ } /*}}}*/ /*FUNCTION PentaRef::GetBprimePattyn {{{*/ -void PentaRef::GetBprimePattyn(double* B, double* xyz_list, GaussPenta* gauss_coord){ +void PentaRef::GetBprimePattyn(IssmDouble* B, IssmDouble* xyz_list, GaussPenta* gauss_coord){ /*Compute B prime matrix. B=[B1 B2 B3 B4 B5 B6] where Bi is of size 5*NDOF2. * For node i, Bi can be expressed in the actual coordinate system * by: @@ -187,7 +187,7 @@ * * We assume B has been allocated already, of size: 5x(NDOF2*NUMNODESP1) */ - double dbasis[3][NUMNODESP1]; + IssmDouble dbasis[3][NUMNODESP1]; /*Get dbasis in actual coordinate system: */ GetNodalFunctionsP1Derivatives(&dbasis[0][0],xyz_list, gauss_coord); @@ -212,7 +212,7 @@ } /*}}}*/ /*FUNCTION PentaRef::GetBprimeMacAyealStokes{{{*/ -void PentaRef::GetBprimeMacAyealStokes(double* Bprime, double* xyz_list, GaussPenta* gauss){ +void PentaRef::GetBprimeMacAyealStokes(IssmDouble* Bprime, IssmDouble* xyz_list, GaussPenta* gauss){ /*Compute Bprime matrix. Bprime=[Bprime1 Bprime2 Bprime3 Bprime4 Bprime5 Bprime6] where Bprimei is of size 5*NDOF2. * For node i, Bprimei can be expressed in the actual coordinate system * by: @@ -225,7 +225,7 @@ */ int i; - double dh1dh7[3][NUMNODESMINI]; + IssmDouble dh1dh7[3][NUMNODESMINI]; /*Get dh1dh6 in actual coordinate system: */ GetNodalFunctionsMINIDerivatives(&dh1dh7[0][0],xyz_list, gauss); @@ -252,7 +252,7 @@ } /*}}}*/ /*FUNCTION PentaRef::GetBStokes {{{*/ -void PentaRef::GetBStokes(double* B, double* xyz_list, GaussPenta* gauss){ +void PentaRef::GetBStokes(IssmDouble* B, IssmDouble* xyz_list, GaussPenta* gauss){ /*Compute B matrix. B=[B1 B2 B3 B4 B5 B6] where Bi is of size 3*NDOF4. * For node i, Bi can be expressed in the actual coordinate system @@ -270,8 +270,8 @@ int i; - double dh1dh7[3][NUMNODESMINI]; - double l1l6[NUMNODESP1]; + IssmDouble dh1dh7[3][NUMNODESMINI]; + IssmDouble l1l6[NUMNODESP1]; /*Get dh1dh7 in actual coordinate system: */ GetNodalFunctionsMINIDerivatives(&dh1dh7[0][0],xyz_list, gauss); @@ -319,7 +319,7 @@ } /*}}}*/ /*FUNCTION PentaRef::GetBprimeStokes {{{*/ -void PentaRef::GetBprimeStokes(double* B_prime, double* xyz_list, GaussPenta* gauss){ +void PentaRef::GetBprimeStokes(IssmDouble* B_prime, IssmDouble* xyz_list, GaussPenta* gauss){ /* Compute B' matrix. B'=[B1' B2' B3' B4' B5' B6' Bb'] where Bi' is of size 3*NDOF2. * For node i, Bi' can be expressed in the actual coordinate system * by: @@ -337,8 +337,8 @@ */ int i; - double dh1dh7[3][NUMNODESMINI]; - double l1l6[NUMNODESP1]; + IssmDouble dh1dh7[3][NUMNODESMINI]; + IssmDouble l1l6[NUMNODESP1]; /*Get dh1dh7 in actual coordinate system: */ GetNodalFunctionsMINIDerivatives(&dh1dh7[0][0],xyz_list, gauss); @@ -386,7 +386,7 @@ } /*}}}*/ /*FUNCTION PentaRef::GetBAdvec{{{*/ -void PentaRef::GetBAdvec(double* B_advec, double* xyz_list, GaussPenta* gauss){ +void PentaRef::GetBAdvec(IssmDouble* B_advec, IssmDouble* xyz_list, GaussPenta* gauss){ /*Compute B matrix. B=[B1 B2 B3 B4 B5 B6] where Bi is of size 5*NDOF1. * For node i, Bi' can be expressed in the actual coordinate system * by: @@ -399,7 +399,7 @@ */ /*Same thing in the actual coordinate system: */ - double l1l6[6]; + IssmDouble l1l6[6]; /*Get dh1dh2dh3 in actual coordinates system : */ GetNodalFunctionsP1(l1l6, gauss); @@ -413,7 +413,7 @@ } /*}}}*/ /*FUNCTION PentaRef::GetBConduct{{{*/ -void PentaRef::GetBConduct(double* B_conduct, double* xyz_list, GaussPenta* gauss){ +void PentaRef::GetBConduct(IssmDouble* B_conduct, IssmDouble* xyz_list, GaussPenta* gauss){ /*Compute B matrix. B=[B1 B2 B3 B4 B5 B6] where Bi is of size 5*NDOF1. * For node i, Bi' can be expressed in the actual coordinate system * by: @@ -426,7 +426,7 @@ */ /*Same thing in the actual coordinate system: */ - double dh1dh6[3][NUMNODESP1]; + IssmDouble dh1dh6[3][NUMNODESP1]; /*Get dh1dh2dh3 in actual coordinates system : */ GetNodalFunctionsP1Derivatives(&dh1dh6[0][0],xyz_list,gauss); @@ -440,12 +440,12 @@ } /*}}}*/ /*FUNCTION PentaRef::GetBVert{{{*/ -void PentaRef::GetBVert(double* B, double* xyz_list, GaussPenta* gauss){ +void PentaRef::GetBVert(IssmDouble* B, IssmDouble* xyz_list, GaussPenta* gauss){ /* Compute B matrix. B=[dh1/dz dh2/dz dh3/dz dh4/dz dh5/dz dh6/dz]; where hi is the interpolation function for node i.*/ int i; - double dh1dh6[3][NUMNODESP1]; + IssmDouble dh1dh6[3][NUMNODESP1]; /*Get dh1dh6 in actual coordinate system: */ GetNodalFunctionsP1Derivatives(&dh1dh6[0][0],xyz_list, gauss); @@ -458,7 +458,7 @@ } /*}}}*/ /*FUNCTION PentaRef::GetBprimeAdvec{{{*/ -void PentaRef::GetBprimeAdvec(double* Bprime_advec, double* xyz_list, GaussPenta* gauss){ +void PentaRef::GetBprimeAdvec(IssmDouble* Bprime_advec, IssmDouble* xyz_list, GaussPenta* gauss){ /*Compute B matrix. B=[B1 B2 B3 B4 B5 B6] where Bi is of size 5*NDOF1. * For node i, Bi' can be expressed in the actual coordinate system * by: @@ -471,7 +471,7 @@ */ /*Same thing in the actual coordinate system: */ - double dh1dh6[3][NUMNODESP1]; + IssmDouble dh1dh6[3][NUMNODESP1]; /*Get dh1dh2dh3 in actual coordinates system : */ GetNodalFunctionsP1Derivatives(&dh1dh6[0][0],xyz_list,gauss); @@ -485,7 +485,7 @@ } /*}}}*/ /*FUNCTION PentaRef::GetBprimeVert{{{*/ -void PentaRef::GetBprimeVert(double* B, double* xyz_list, GaussPenta* gauss){ +void PentaRef::GetBprimeVert(IssmDouble* B, IssmDouble* xyz_list, GaussPenta* gauss){ /* Compute Bprime matrix. Bprime=[L1 L2 L3 L4 L5 L6] where Li is the nodal function for node i*/ GetNodalFunctionsP1(B, gauss); @@ -493,7 +493,7 @@ } /*}}}*/ /*FUNCTION PentaRef::GetL{{{*/ -void PentaRef::GetL(double* L, GaussPenta* gauss, int numdof){ +void PentaRef::GetL(IssmDouble* L, GaussPenta* gauss, int numdof){ /*Compute L matrix. L=[L1 L2 L3] where Li is square and of size numdof. ** For node i, Li can be expressed in the actual coordinate system ** by: @@ -508,7 +508,7 @@ **/ int i; - double l1l6[6]; + IssmDouble l1l6[6]; /*Get l1l6 in actual coordinate system: */ GetNodalFunctionsP1(l1l6,gauss); @@ -530,7 +530,7 @@ } /*}}}*/ /*FUNCTION PentaRef::GetLStokes{{{*/ -void PentaRef::GetLStokes(double* LStokes, GaussPenta* gauss){ +void PentaRef::GetLStokes(IssmDouble* LStokes, GaussPenta* gauss){ /* * Compute L matrix. L=[L1 L2 L3] where Li is square and of size numdof. * For node i, Li can be expressed in the actual coordinate system @@ -543,7 +543,7 @@ */ const int num_dof=4; - double l1l2l3[NUMNODESP1_2d]; + IssmDouble l1l2l3[NUMNODESP1_2d]; /*Get l1l2l3 in actual coordinate system: */ l1l2l3[0]=gauss->coord1*(1-gauss->coord4)/2.0; @@ -565,7 +565,7 @@ } /*}}}*/ /*FUNCTION PentaRef::GetLprimeStokes {{{*/ -void PentaRef::GetLprimeStokes(double* LprimeStokes, double* xyz_list, GaussPenta* gauss){ +void PentaRef::GetLprimeStokes(IssmDouble* LprimeStokes, IssmDouble* xyz_list, GaussPenta* gauss){ /* * Compute Lprime matrix. Lprime=[Lp1 Lp2 Lp3] where Lpi is square and of size numdof. @@ -605,8 +605,8 @@ int i; int num_dof=4; - double l1l2l3[NUMNODESP1_2d]; - double dh1dh6[3][NUMNODESP1]; + IssmDouble l1l2l3[NUMNODESP1_2d]; + IssmDouble dh1dh6[3][NUMNODESP1]; /*Get l1l2l3 in actual coordinate system: */ l1l2l3[0]=gauss->coord1*(1-gauss->coord4)/2.0; @@ -677,7 +677,7 @@ } /*}}}*/ /*FUNCTION PentaRef::GetLMacAyealStokes {{{*/ -void PentaRef::GetLMacAyealStokes(double* LStokes, GaussPenta* gauss){ +void PentaRef::GetLMacAyealStokes(IssmDouble* LStokes, GaussPenta* gauss){ /* * Compute L matrix. L=[L1 L2 L3] where Li is square and of size numdof. * For node i, Li can be expressed in the actual coordinate system @@ -696,7 +696,7 @@ int i; int num_dof=2; - double l1l2l3[NUMNODESP1_2d]; + IssmDouble l1l2l3[NUMNODESP1_2d]; /*Get l1l2l3 in actual coordinate system: */ @@ -727,7 +727,7 @@ } /*}}}*/ /*FUNCTION PentaRef::GetLprimeMacAyealStokes {{{*/ -void PentaRef::GetLprimeMacAyealStokes(double* LprimeStokes, double* xyz_list, GaussPenta* gauss){ +void PentaRef::GetLprimeMacAyealStokes(IssmDouble* LprimeStokes, IssmDouble* xyz_list, GaussPenta* gauss){ /* * Compute Lprime matrix. Lprime=[Lp1 Lp2 Lp3] where Lpi is square and of size numdof. @@ -746,8 +746,8 @@ int i; int num_dof=4; - double l1l2l3[NUMNODESP1_2d]; - double dh1dh6[3][NUMNODESP1]; + IssmDouble l1l2l3[NUMNODESP1_2d]; + IssmDouble dh1dh6[3][NUMNODESP1]; /*Get l1l2l3 in actual coordinate system: */ l1l2l3[0]=gauss->coord1*(1-gauss->coord4)/2.0; @@ -794,7 +794,7 @@ } /*}}}*/ /*FUNCTION PentaRef::GetLStokesMacAyeal {{{*/ -void PentaRef::GetLStokesMacAyeal(double* LStokes, GaussPenta* gauss){ +void PentaRef::GetLStokesMacAyeal(IssmDouble* LStokes, GaussPenta* gauss){ /* * Compute L matrix. L=[L1 L2 L3] where Li is square and of size numdof. * For node i, Li can be expressed in the actual coordinate system @@ -809,7 +809,7 @@ int i; int num_dof=4; - double l1l2l3[NUMNODESP1_2d]; + IssmDouble l1l2l3[NUMNODESP1_2d]; /*Get l1l2l3 in actual coordinate system: */ @@ -840,7 +840,7 @@ } /*}}}*/ /*FUNCTION PentaRef::GetLprimeStokesMacAyeal {{{*/ -void PentaRef::GetLprimeStokesMacAyeal(double* LprimeStokes, double* xyz_list, GaussPenta* gauss){ +void PentaRef::GetLprimeStokesMacAyeal(IssmDouble* LprimeStokes, IssmDouble* xyz_list, GaussPenta* gauss){ /* * Compute Lprime matrix. Lprime=[Lp1 Lp2 Lp3] where Lpi is square and of size numdof. @@ -855,8 +855,8 @@ int i; int num_dof=2; - double l1l2l3[NUMNODESP1_2d]; - double dh1dh6[3][NUMNODESP1]; + IssmDouble l1l2l3[NUMNODESP1_2d]; + IssmDouble dh1dh6[3][NUMNODESP1]; /*Get l1l2l3 in actual coordinate system: */ l1l2l3[0]=gauss->coord1*(1-gauss->coord4)/2.0; @@ -879,19 +879,19 @@ } /*}}}*/ /*FUNCTION PentaRef::GetJacobian {{{*/ -void PentaRef::GetJacobian(double* J, double* xyz_list,GaussPenta* gauss){ +void PentaRef::GetJacobian(IssmDouble* J, IssmDouble* xyz_list,GaussPenta* gauss){ int i,j; /*The Jacobian is constant over the element, discard the gaussian points. * J is assumed to have been allocated of size NDOF2xNDOF2.*/ - double A1,A2,A3; //area coordinates - double xi,eta,zi; //parametric coordinates + IssmDouble A1,A2,A3; //area coordinates + IssmDouble xi,eta,zi; //parametric coordinates - double x1,x2,x3,x4,x5,x6; - double y1,y2,y3,y4,y5,y6; - double z1,z2,z3,z4,z5,z6; + IssmDouble x1,x2,x3,x4,x5,x6; + IssmDouble y1,y2,y3,y4,y5,y6; + IssmDouble z1,z2,z3,z4,z5,z6; /*Figure out xi,eta and zi (parametric coordinates), for this gaussian point: */ A1=gauss->coord1; @@ -938,10 +938,10 @@ } /*}}}*/ /*FUNCTION PentaRef::GetJacobianDeterminant {{{*/ -void PentaRef::GetJacobianDeterminant(double* Jdet, double* xyz_list,GaussPenta* gauss){ +void PentaRef::GetJacobianDeterminant(IssmDouble* Jdet, IssmDouble* xyz_list,GaussPenta* gauss){ /*On a penta, Jacobian varies according to coordinates. We need to get the Jacobian, and take * the determinant of it: */ - double J[3][3]; + IssmDouble J[3][3]; /*Get Jacobian*/ GetJacobian(&J[0][0],xyz_list,gauss); @@ -953,11 +953,11 @@ } /*}}}*/ /*FUNCTION PentaRef::GetTriaJacobianDeterminant{{{*/ -void PentaRef::GetTriaJacobianDeterminant(double* Jdet, double* xyz_list,GaussPenta* gauss){ +void PentaRef::GetTriaJacobianDeterminant(IssmDouble* Jdet, IssmDouble* xyz_list,GaussPenta* gauss){ /*The Jacobian determinant is constant over the element, discard the gaussian points. * J is assumed to have been allocated of size NDOF2xNDOF2.*/ - double x1,x2,x3,y1,y2,y3,z1,z2,z3; + IssmDouble x1,x2,x3,y1,y2,y3,z1,z2,z3; x1=*(xyz_list+3*0+0); y1=*(xyz_list+3*0+1); @@ -975,11 +975,11 @@ } /*}}}*/ /*FUNCTION PentaRef::GetSegmentJacobianDeterminant{{{*/ -void PentaRef::GetSegmentJacobianDeterminant(double* Jdet, double* xyz_list,GaussPenta* gauss){ +void PentaRef::GetSegmentJacobianDeterminant(IssmDouble* Jdet, IssmDouble* xyz_list,GaussPenta* gauss){ /*The Jacobian determinant is constant over the element, discard the gaussian points. * J is assumed to have been allocated of size NDOF2xNDOF2.*/ - double x1,x2,y1,y2,z1,z2; + IssmDouble x1,x2,y1,y2,z1,z2; x1=*(xyz_list+3*0+0); y1=*(xyz_list+3*0+1); @@ -994,10 +994,10 @@ } /*}}}*/ /*FUNCTION PentaRef::GetJacobianInvert {{{*/ -void PentaRef::GetJacobianInvert(double* Jinv, double* xyz_list,GaussPenta* gauss){ +void PentaRef::GetJacobianInvert(IssmDouble* Jinv, IssmDouble* xyz_list,GaussPenta* gauss){ /*Jacobian*/ - double J[3][3]; + IssmDouble J[3][3]; /*Call Jacobian routine to get the jacobian:*/ GetJacobian(&J[0][0], xyz_list, gauss); @@ -1007,7 +1007,7 @@ } /*}}}*/ /*FUNCTION PentaRef::GetNodalFunctionsMINI{{{*/ -void PentaRef::GetNodalFunctionsMINI(double* l1l7, GaussPenta* gauss){ +void PentaRef::GetNodalFunctionsMINI(IssmDouble* l1l7, GaussPenta* gauss){ /*This routine returns the values of the nodal functions at the gaussian point.*/ l1l7[0]=gauss->coord1*(1.0-gauss->coord4)/2.0; @@ -1021,14 +1021,14 @@ } /*}}}*/ /*FUNCTION PentaRef::GetNodalFunctionsMINIDerivatives{{{*/ -void PentaRef::GetNodalFunctionsMINIDerivatives(double* dh1dh7,double* xyz_list, GaussPenta* gauss){ +void PentaRef::GetNodalFunctionsMINIDerivatives(IssmDouble* dh1dh7,IssmDouble* xyz_list, GaussPenta* gauss){ /*This routine returns the values of the nodal functions derivatives (with respect to the * actual coordinate system): */ int i; - double dh1dh7_ref[3][NUMNODESMINI]; - double Jinv[3][3]; + IssmDouble dh1dh7_ref[3][NUMNODESMINI]; + IssmDouble Jinv[3][3]; /*Get derivative values with respect to parametric coordinate system: */ GetNodalFunctionsMINIDerivativesReference(&dh1dh7_ref[0][0], gauss); @@ -1052,13 +1052,13 @@ } /*}}}*/ /*FUNCTION PentaRef::GetNodalFunctionsMINIDerivativesReference{{{*/ -void PentaRef::GetNodalFunctionsMINIDerivativesReference(double* dl1dl7,GaussPenta* gauss){ +void PentaRef::GetNodalFunctionsMINIDerivativesReference(IssmDouble* dl1dl7,GaussPenta* gauss){ /*This routine returns the values of the nodal functions derivatives (with respect to the * natural coordinate system) at the gaussian point. */ - double r=gauss->coord2-gauss->coord1; - double s=-3.0/SQRT3*(gauss->coord1+gauss->coord2-2.0/3.0); - double zeta=gauss->coord4; + IssmDouble r=gauss->coord2-gauss->coord1; + IssmDouble s=-3.0/SQRT3*(gauss->coord1+gauss->coord2-2.0/3.0); + IssmDouble zeta=gauss->coord4; /*First nodal function: */ *(dl1dl7+NUMNODESMINI*0+0)=-0.5*(1.0-zeta)/2.0; @@ -1098,7 +1098,7 @@ } /*}}}*/ /*FUNCTION PentaRef::GetNodalFunctionsP1 {{{*/ -void PentaRef::GetNodalFunctionsP1(double* l1l6, GaussPenta* gauss){ +void PentaRef::GetNodalFunctionsP1(IssmDouble* l1l6, GaussPenta* gauss){ /*This routine returns the values of the nodal functions at the gaussian point.*/ l1l6[0]=gauss->coord1*(1-gauss->coord4)/2.0; @@ -1111,12 +1111,12 @@ } /*}}}*/ /*FUNCTION PentaRef::GetNodalFunctionsP1Derivatives {{{*/ -void PentaRef::GetNodalFunctionsP1Derivatives(double* dh1dh6,double* xyz_list, GaussPenta* gauss){ +void PentaRef::GetNodalFunctionsP1Derivatives(IssmDouble* dh1dh6,IssmDouble* xyz_list, GaussPenta* gauss){ /*This routine returns the values of the nodal functions derivatives (with respect to the * actual coordinate system): */ - double dh1dh6_ref[NDOF3][NUMNODESP1]; - double Jinv[NDOF3][NDOF3]; + IssmDouble dh1dh6_ref[NDOF3][NUMNODESP1]; + IssmDouble Jinv[NDOF3][NDOF3]; /*Get derivative values with respect to parametric coordinate system: */ GetNodalFunctionsP1DerivativesReference(&dh1dh6_ref[0][0], gauss); @@ -1140,12 +1140,12 @@ } /*}}}*/ /*FUNCTION PentaRef::GetNodalFunctionsP1DerivativesReference {{{*/ -void PentaRef::GetNodalFunctionsP1DerivativesReference(double* dl1dl6,GaussPenta* gauss){ +void PentaRef::GetNodalFunctionsP1DerivativesReference(IssmDouble* dl1dl6,GaussPenta* gauss){ /*This routine returns the values of the nodal functions derivatives (with respect to the * natural coordinate system) at the gaussian point. Those values vary along xi,eta,z */ - double A1,A2,A3,z; + IssmDouble A1,A2,A3,z; A1=gauss->coord1; _assert_(A1>=0 && A1<=1);//first area coordinate value. In term of xi and eta: A1=(1-xi)/2-eta/(2*SQRT3); A2=gauss->coord2; _assert_(A2>=0 && A2<=1);//second area coordinate value In term of xi and eta: A2=(1+xi)/2-eta/(2*SQRT3); @@ -1184,10 +1184,10 @@ } /*}}}*/ /*FUNCTION PentaRef::GetQuadNodalFunctions {{{*/ -void PentaRef::GetQuadNodalFunctions(double* l1l4,GaussPenta* gauss,int index1,int index2,int index3,int index4){ +void PentaRef::GetQuadNodalFunctions(IssmDouble* l1l4,GaussPenta* gauss,int index1,int index2,int index3,int index4){ /*This routine returns the values of the nodal functions at the gaussian point.*/ - double BasisFunctions[6]; + IssmDouble BasisFunctions[6]; GetNodalFunctionsP1(&BasisFunctions[0],gauss); @@ -1204,10 +1204,10 @@ } /*}}}*/ /*FUNCTION PentaRef::GetQuadJacobianDeterminant{{{*/ -void PentaRef::GetQuadJacobianDeterminant(double* Jdet,double xyz_list[4][3],GaussPenta* gauss){ +void PentaRef::GetQuadJacobianDeterminant(IssmDouble* Jdet,IssmDouble xyz_list[4][3],GaussPenta* gauss){ /*This routine returns the values of the nodal functions at the gaussian point.*/ - double x1,x2,x3,x4,y1,y2,y3,y4,z1,z2,z3,z4; + IssmDouble x1,x2,x3,x4,y1,y2,y3,y4,z1,z2,z3,z4; x1=xyz_list[0][0]; y1=xyz_list[0][1]; @@ -1230,11 +1230,11 @@ } /*}}}*/ /*FUNCTION PentaRef::GetInputValue{{{*/ -void PentaRef::GetInputValue(double* pvalue,double* plist,GaussPenta* gauss){ +void PentaRef::GetInputValue(IssmDouble* pvalue,IssmDouble* plist,GaussPenta* gauss){ /*P1 interpolation on Gauss point*/ /*intermediary*/ - double l1l6[6]; + IssmDouble l1l6[6]; /*nodal functions: */ GetNodalFunctionsP1(&l1l6[0],gauss); @@ -1245,7 +1245,7 @@ } /*}}}*/ /*FUNCTION PentaRef::GetInputDerivativeValue{{{*/ -void PentaRef::GetInputDerivativeValue(double* p, double* plist,double* xyz_list, GaussPenta* gauss){ +void PentaRef::GetInputDerivativeValue(IssmDouble* p, IssmDouble* plist,IssmDouble* xyz_list, GaussPenta* gauss){ /*From node values of parameter p (p_list[0], p_list[1], p_list[2], p_list[3], p_list[4] and p_list[4]), return parameter derivative value at gaussian point specified by gauss_coord: * dp/dx=p_list[0]*dh1/dx+p_list[1]*dh2/dx+p_list[2]*dh3/dx+p_list[3]*dh4/dx+p_list[4]*dh5/dx+p_list[5]*dh6/dx; * dp/dy=p_list[0]*dh1/dy+p_list[1]*dh2/dy+p_list[2]*dh3/dy+p_list[3]*dh4/dy+p_list[4]*dh5/dy+p_list[5]*dh6/dy; @@ -1253,7 +1253,7 @@ * * p is a vector of size 3x1 already allocated. */ - double dh1dh6[3][NUMNODESP1]; + IssmDouble dh1dh6[3][NUMNODESP1]; /*Get nodal funnctions derivatives in actual coordinate system: */ GetNodalFunctionsP1Derivatives(&dh1dh6[0][0],xyz_list, gauss); Index: /u/astrid-r1b/morlighe/issmuci/trunk-jpl/../trunk-jpl/src/c/objects/Elements/Tria.h =================================================================== --- /u/astrid-r1b/morlighe/issmuci/trunk-jpl/../trunk-jpl/src/c/objects/Elements/Tria.h (revision 12470) +++ /u/astrid-r1b/morlighe/issmuci/trunk-jpl/../trunk-jpl/src/c/objects/Elements/Tria.h (revision 12471) @@ -54,23 +54,23 @@ Object* copy(); /*}}}*/ /*Update virtual functions resolution: {{{*/ - void InputUpdateFromSolution(double* solutiong); - void InputUpdateFromVector(double* vector, int name, int type); + void InputUpdateFromSolution(IssmDouble* solutiong); + void InputUpdateFromVector(IssmDouble* vector, int name, int type); void InputUpdateFromVector(int* vector, int name, int type); void InputUpdateFromVector(bool* vector, int name, int type); #ifdef _HAVE_DAKOTA_ - void InputUpdateFromVectorDakota(double* vector, int name, int type); + void InputUpdateFromVectorDakota(IssmDouble* vector, int name, int type); void InputUpdateFromVectorDakota(int* vector, int name, int type); void InputUpdateFromVectorDakota(bool* vector, int name, int type); - void InputUpdateFromMatrixDakota(double* matrix, int nows, int ncols, int name, int type); + void InputUpdateFromMatrixDakota(IssmDouble* matrix, int nows, int ncols, int name, int type); #endif - void InputUpdateFromConstant(double constant, int name); + void InputUpdateFromConstant(IssmDouble constant, int name); void InputUpdateFromConstant(int constant, int name); void InputUpdateFromConstant(bool constant, int name); void InputUpdateFromIoModel(int index, IoModel* iomodel); /*}}}*/ /*Element virtual functions definitions: {{{*/ - void AverageOntoPartition(Vector* partition_contributions,Vector* partition_areas,double* vertex_response,double* qmu_part); + void AverageOntoPartition(Vector* partition_contributions,Vector* partition_areas,IssmDouble* vertex_response,IssmDouble* qmu_part); void ComputeBasalStress(Vector* sigma_b); void ComputeStrainRate(Vector* eps); void ComputeStressTensor(); @@ -84,58 +84,58 @@ bool IsOnBed(); bool IsFloating(); bool IsNodeOnShelf(); - bool IsNodeOnShelfFromFlags(double* flags); + bool IsNodeOnShelfFromFlags(IssmDouble* flags); bool IsOnWater(); void GetSolutionFromInputs(Vector* solution); void GetVectorFromInputs(Vector* vector, int name_enum); void GetVectorFromResults(Vector* vector,int offset,int interp); - void InputArtificialNoise(int enum_type,double min, double max); - bool InputConvergence(double* eps, int* enums,int num_enums,int* criterionenums,double* criterionvalues,int num_criterionenums); - void InputCreate(double scalar,int name,int code); - void InputCreate(double* vector, int index,IoModel* iomodel,int M,int N,int vector_type,int vector_enum,int code); + void InputArtificialNoise(int enum_type,IssmDouble min, IssmDouble max); + bool InputConvergence(IssmDouble* eps, int* enums,int num_enums,int* criterionenums,IssmDouble* criterionvalues,int num_criterionenums); + void InputCreate(IssmDouble scalar,int name,int code); + void InputCreate(IssmDouble* vector, int index,IoModel* iomodel,int M,int N,int vector_type,int vector_enum,int code); void InputDepthAverageAtBase(int enum_type,int average_enum_type,int object_enum=MeshElementsEnum); void InputDuplicate(int original_enum,int new_enum); - void InputScale(int enum_type,double scale_factor); - void InputToResult(int enum_type,int step,double time); + void InputScale(int enum_type,IssmDouble scale_factor); + void InputToResult(int enum_type,int step,IssmDouble time); void DeleteResults(void); void MaterialUpdateFromTemperature(void){_error_("not implemented yet");}; - void MigrateGroundingLine(double* oldfloating,double* sheet_ungrounding); - int NodalValue(double* pvalue, int index, int natureofdataenum,bool process_units); + void MigrateGroundingLine(IssmDouble* oldfloating,IssmDouble* sheet_ungrounding); + int NodalValue(IssmDouble* pvalue, int index, int natureofdataenum,bool process_units); void PotentialSheetUngrounding(Vector* potential_sheet_ungrounding); - void PositiveDegreeDay(double* pdds,double* pds,double signorm); - void RequestedOutput(int output_enum,int step,double time); - void ListResultsInfo(int** results_enums,int** results_size,double** results_times,int** results_steps,int* num_results); + void PositiveDegreeDay(IssmDouble* pdds,IssmDouble* pds,IssmDouble signorm); + void RequestedOutput(int output_enum,int step,IssmDouble time); + void ListResultsInfo(int** results_enums,int** results_size,IssmDouble** results_times,int** results_steps,int* num_results); void PatchFill(int* pcount, Patch* patch); void PatchSize(int* pnumrows, int* pnumvertices,int* pnumnodes); void ProcessResultsUnits(void); void ResetCoordinateSystem(void){_error_("not implemented yet");}; - double SurfaceArea(void); + IssmDouble SurfaceArea(void); void Update(int index, IoModel* iomodel,int analysis_counter,int analysis_type); - int UpdatePotentialSheetUngrounding(double* vertices_potentially_ungrounding,Vector* vec_nodes_on_iceshelf,double* nodes_on_iceshelf); - double TimeAdapt(); + int UpdatePotentialSheetUngrounding(IssmDouble* vertices_potentially_ungrounding,Vector* vec_nodes_on_iceshelf,IssmDouble* nodes_on_iceshelf); + IssmDouble TimeAdapt(); int* GetHorizontalNeighboorSids(void); - void SmearFunction(Vector* smearedvector,double (*WeightFunction)(double distance,double radius),double radius); + void SmearFunction(Vector* smearedvector,IssmDouble (*WeightFunction)(IssmDouble distance,IssmDouble radius),IssmDouble radius); #ifdef _HAVE_RESPONSES_ - double IceVolume(void); - 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 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 ElementResponse(double* presponse,int response_enum,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); + IssmDouble IceVolume(void); + void MinVel(IssmDouble* pminvel, bool process_units); + void MinVx(IssmDouble* pminvx, bool process_units); + void MinVy(IssmDouble* pminvy, bool process_units); + void MinVz(IssmDouble* pminvz, bool process_units); + IssmDouble MassFlux(IssmDouble* segment,bool process_units); + void MaxAbsVx(IssmDouble* pmaxabsvx, bool process_units); + void MaxAbsVy(IssmDouble* pmaxabsvy, bool process_units); + void MaxAbsVz(IssmDouble* pmaxabsvz, bool process_units); + void ElementResponse(IssmDouble* presponse,int response_enum,bool process_units); + void MaxVel(IssmDouble* pmaxvel, bool process_units); + void MaxVx(IssmDouble* pmaxvx, bool process_units); + void MaxVy(IssmDouble* pmaxvy, bool process_units); + void MaxVz(IssmDouble* pmaxvz, bool process_units); #endif #ifdef _HAVE_CONTROL_ - double DragCoefficientAbsGradient(bool process_units,int weight_index); + IssmDouble DragCoefficientAbsGradient(bool process_units,int weight_index); void GradientIndexing(int* indexing,int control_index); void Gradj(Vector* gradient,int control_type,int control_index); void GradjBGradient(Vector* gradient,int weight_index,int control_index); @@ -147,19 +147,19 @@ void GradjVxBalancedthickness(Vector* gradient,int control_index); void GradjVyBalancedthickness(Vector* gradient,int control_index); void GetVectorFromControlInputs(Vector* gradient,int control_enum,int control_index,const char* data); - void SetControlInputsFromVector(double* vector,int control_enum,int control_index); + void SetControlInputsFromVector(IssmDouble* vector,int control_enum,int control_index); void ControlInputGetGradient(Vector* gradient,int enum_type,int control_index); - void ControlInputScaleGradient(int enum_type,double scale); - void ControlInputSetGradient(double* gradient,int enum_type,int control_index); - double RheologyBbarAbsGradient(bool process_units,int weight_index); - double ThicknessAbsMisfit( bool process_units,int weight_index); - double SurfaceAbsVelMisfit( bool process_units,int weight_index); - double ThicknessAbsGradient(bool process_units,int weight_index); - double SurfaceRelVelMisfit( bool process_units,int weight_index); - double SurfaceLogVelMisfit( bool process_units,int weight_index); - double SurfaceLogVxVyMisfit( bool process_units,int weight_index); - double SurfaceAverageVelMisfit(bool process_units,int weight_index); - void InputControlUpdate(double scalar,bool save_parameter); + void ControlInputScaleGradient(int enum_type,IssmDouble scale); + void ControlInputSetGradient(IssmDouble* gradient,int enum_type,int control_index); + IssmDouble RheologyBbarAbsGradient(bool process_units,int weight_index); + IssmDouble ThicknessAbsMisfit( bool process_units,int weight_index); + IssmDouble SurfaceAbsVelMisfit( bool process_units,int weight_index); + IssmDouble ThicknessAbsGradient(bool process_units,int weight_index); + IssmDouble SurfaceRelVelMisfit( bool process_units,int weight_index); + IssmDouble SurfaceLogVelMisfit( bool process_units,int weight_index); + IssmDouble SurfaceLogVxVyMisfit( bool process_units,int weight_index); + IssmDouble SurfaceAverageVelMisfit(bool process_units,int weight_index); + void InputControlUpdate(IssmDouble scalar,bool save_parameter); #endif /*}}}*/ @@ -179,22 +179,22 @@ ElementVector* CreatePVectorPrognostic_CG(void); ElementVector* CreatePVectorPrognostic_DG(void); ElementVector* CreatePVectorSlope(void); - double GetArea(void); + IssmDouble GetArea(void); int GetElementType(void); void GetDofList(int** pdoflist,int approximation_enum,int setenum); void GetDofList1(int* doflist); void GetSidList(int* sidlist); void GetConnectivityList(int* connectivity); - void GetInputListOnVertices(double* pvalue,int enumtype); - void GetInputListOnVertices(double* pvalue,int enumtype,double defaultvalue); - void GetInputListOnVertices(double* pvalue,int enumtype,double defaultvalue,int index); //TO BE REMOVED - void GetInputValue(double* pvalue,Node* node,int enumtype); - void GetStrainRate2d(double* epsilon,double* xyz_list, GaussTria* gauss, Input* vx_input, Input* vy_input); - void InputUpdateFromSolutionOneDof(double* solution,int enum_type); - void InputUpdateFromSolutionPrognostic(double* solution); + void GetInputListOnVertices(IssmDouble* pvalue,int enumtype); + void GetInputListOnVertices(IssmDouble* pvalue,int enumtype,IssmDouble defaultvalue); + void GetInputListOnVertices(IssmDouble* pvalue,int enumtype,IssmDouble defaultvalue,int index); //TO BE REMOVED + void GetInputValue(IssmDouble* pvalue,Node* node,int enumtype); + void GetStrainRate2d(IssmDouble* epsilon,IssmDouble* xyz_list, GaussTria* gauss, Input* vx_input, Input* vy_input); + void InputUpdateFromSolutionOneDof(IssmDouble* solution,int enum_type); + void InputUpdateFromSolutionPrognostic(IssmDouble* solution); bool IsInput(int name); void SetClone(int* minranks); - void SurfaceNormal(double* surface_normal, double xyz_list[3][3]); + void SurfaceNormal(IssmDouble* surface_normal, IssmDouble xyz_list[3][3]); #ifdef _HAVE_DIAGNOSTIC_ ElementMatrix* CreateKMatrixDiagnosticMacAyeal(void); @@ -206,8 +206,8 @@ ElementMatrix* CreateJacobianDiagnosticMacayeal(void); void GetSolutionFromInputsDiagnosticHoriz(Vector* solution); void GetSolutionFromInputsDiagnosticHutter(Vector* solution); - void InputUpdateFromSolutionDiagnosticHoriz( double* solution); - void InputUpdateFromSolutionDiagnosticHutter( double* solution); + void InputUpdateFromSolutionDiagnosticHoriz( IssmDouble* solution); + void InputUpdateFromSolutionDiagnosticHutter( IssmDouble* solution); #endif #ifdef _HAVE_CONTROL_ @@ -216,8 +216,8 @@ ElementVector* CreatePVectorAdjointHoriz(void); ElementVector* CreatePVectorAdjointStokes(void); ElementVector* CreatePVectorAdjointBalancethickness(void); - void InputUpdateFromSolutionAdjointBalancethickness( double* solution); - void InputUpdateFromSolutionAdjointHoriz( double* solution); + void InputUpdateFromSolutionAdjointBalancethickness( IssmDouble* solution); + void InputUpdateFromSolutionAdjointHoriz( IssmDouble* solution); #endif #ifdef _HAVE_HYDROLOGY_ @@ -225,7 +225,7 @@ ElementVector* CreatePVectorHydrology(void); void CreateHydrologyWaterVelocityInput(void); void GetSolutionFromInputsHydrology(Vector* solution); - void InputUpdateFromSolutionHydrology(double* solution); + void InputUpdateFromSolutionHydrology(IssmDouble* solution); #endif #ifdef _HAVE_BALANCED_ #endif Index: /u/astrid-r1b/morlighe/issmuci/trunk-jpl/../trunk-jpl/src/c/objects/Elements/Penta.cpp =================================================================== --- /u/astrid-r1b/morlighe/issmuci/trunk-jpl/../trunk-jpl/src/c/objects/Elements/Penta.cpp (revision 12470) +++ /u/astrid-r1b/morlighe/issmuci/trunk-jpl/../trunk-jpl/src/c/objects/Elements/Penta.cpp (revision 12471) @@ -144,17 +144,17 @@ /*Other*/ /*FUNCTION Penta::AverageOntoPartition {{{*/ -void Penta::AverageOntoPartition(Vector* partition_contributions,Vector* partition_areas,double* vertex_response,double* qmu_part){ +void Penta::AverageOntoPartition(Vector* partition_contributions,Vector* partition_areas,IssmDouble* vertex_response,IssmDouble* qmu_part){ _error_("Not supported yet!"); } /*}}}*/ /*FUNCTION Penta::BedNormal {{{*/ -void Penta::BedNormal(double* bed_normal, double xyz_list[3][3]){ +void Penta::BedNormal(IssmDouble* bed_normal, IssmDouble xyz_list[3][3]){ int i; - double v13[3],v23[3]; - double normal[3]; - double normal_norm; + IssmDouble v13[3],v23[3]; + IssmDouble normal[3]; + IssmDouble normal_norm; for (i=0;i<3;i++){ v13[i]=xyz_list[0][i]-xyz_list[2][i]; @@ -180,8 +180,8 @@ /*Intermediaries */ int count,ig; - double basalfriction[NUMVERTICES]={0,0,0,0,0,0}; - double alpha2,vx,vy; + IssmDouble basalfriction[NUMVERTICES]={0,0,0,0,0,0}; + IssmDouble alpha2,vx,vy; Friction* friction=NULL; GaussPenta* gauss=NULL; @@ -232,18 +232,18 @@ int dofp[1]={3}; int analysis_type,approximation; int doflist[NUMVERTICES]; - double xyz_list[NUMVERTICES][3]; - double xyz_list_tria[3][3]; - double rho_ice,gravity,stokesreconditioning; - double pressure,viscosity,bed,Jdet2d; - double bed_normal[3]; - double basalforce[3]; - double epsilon[6]; /* epsilon=[exx,eyy,ezz,exy,exz,eyz];*/ - double devstresstensor[6]; /* epsilon=[exx,eyy,ezz,exy,exz,eyz];*/ - double stresstensor[6]={0.0}; - double sigma_xx,sigma_yy,sigma_zz; - double sigma_xy,sigma_xz,sigma_yz; - double surface=0,value=0; + IssmDouble xyz_list[NUMVERTICES][3]; + IssmDouble xyz_list_tria[3][3]; + IssmDouble rho_ice,gravity,stokesreconditioning; + IssmDouble pressure,viscosity,bed,Jdet2d; + IssmDouble bed_normal[3]; + IssmDouble basalforce[3]; + IssmDouble epsilon[6]; /* epsilon=[exx,eyy,ezz,exy,exz,eyz];*/ + IssmDouble devstresstensor[6]; /* epsilon=[exx,eyy,ezz,exy,exz,eyz];*/ + IssmDouble stresstensor[6]={0.0}; + IssmDouble sigma_xx,sigma_yy,sigma_zz; + IssmDouble sigma_xy,sigma_xz,sigma_yz; + IssmDouble surface=0,value=0; GaussPenta* gauss; /*retrive parameters: */ @@ -325,15 +325,15 @@ void Penta::ComputeStressTensor(){ int iv; - double xyz_list[NUMVERTICES][3]; - double pressure,viscosity; - double epsilon[6]; /* epsilon=[exx,eyy,exy];*/ - double sigma_xx[NUMVERTICES]; - double sigma_yy[NUMVERTICES]; - double sigma_zz[NUMVERTICES]; - double sigma_xy[NUMVERTICES]; - double sigma_xz[NUMVERTICES]; - double sigma_yz[NUMVERTICES]; + IssmDouble xyz_list[NUMVERTICES][3]; + IssmDouble pressure,viscosity; + IssmDouble epsilon[6]; /* epsilon=[exx,eyy,exy];*/ + IssmDouble sigma_xx[NUMVERTICES]; + IssmDouble sigma_yy[NUMVERTICES]; + IssmDouble sigma_zz[NUMVERTICES]; + IssmDouble sigma_xy[NUMVERTICES]; + IssmDouble sigma_xz[NUMVERTICES]; + IssmDouble sigma_yz[NUMVERTICES]; GaussPenta* gauss=NULL; /* Get node coordinates and dof list: */ @@ -763,11 +763,11 @@ } /*}}}*/ /*FUNCTION Penta::GetElementSizes{{{*/ -void Penta::GetElementSizes(double* hx,double* hy,double* hz){ +void Penta::GetElementSizes(IssmDouble* hx,IssmDouble* hy,IssmDouble* hz){ - double xyz_list[NUMVERTICES][3]; - double xmin,ymin,zmin; - double xmax,ymax,zmax; + IssmDouble xyz_list[NUMVERTICES][3]; + IssmDouble xmin,ymin,zmin; + IssmDouble xmax,ymax,zmax; /*Get xyz list: */ GetVerticesCoordinates(&xyz_list[0][0],nodes,NUMVERTICES); @@ -819,11 +819,11 @@ } /*}}}*/ -/*FUNCTION Penta::GetInputListOnVertices(double* pvalue,int enumtype) {{{*/ -void Penta::GetInputListOnVertices(double* pvalue,int enumtype){ +/*FUNCTION Penta::GetInputListOnVertices(IssmDouble* pvalue,int enumtype) {{{*/ +void Penta::GetInputListOnVertices(IssmDouble* pvalue,int enumtype){ /*Intermediaries*/ - double value[NUMVERTICES]; + IssmDouble value[NUMVERTICES]; GaussPenta *gauss = NULL; /*Recover input*/ @@ -844,11 +844,11 @@ delete gauss; } /*}}}*/ -/*FUNCTION Penta::GetInputListOnVertices(double* pvalue,int enumtype,double defaultvalue) {{{*/ -void Penta::GetInputListOnVertices(double* pvalue,int enumtype,double defaultvalue){ +/*FUNCTION Penta::GetInputListOnVertices(IssmDouble* pvalue,int enumtype,IssmDouble defaultvalue) {{{*/ +void Penta::GetInputListOnVertices(IssmDouble* pvalue,int enumtype,IssmDouble defaultvalue){ /*Intermediaries*/ - double value[NUMVERTICES]; + IssmDouble value[NUMVERTICES]; GaussPenta *gauss = NULL; /*Recover input*/ @@ -873,8 +873,8 @@ delete gauss; } /*}}}*/ -/*FUNCTION Penta::GetInputValue(double* pvalue,Node* node,int enumtype) {{{*/ -void Penta::GetInputValue(double* pvalue,Node* node,int enumtype){ +/*FUNCTION Penta::GetInputValue(IssmDouble* pvalue,Node* node,int enumtype) {{{*/ +void Penta::GetInputValue(IssmDouble* pvalue,Node* node,int enumtype){ Input* input=inputs->GetInput(enumtype); if(!input) _error_("No input of type %s found in tria",EnumToStringx(enumtype)); @@ -887,12 +887,12 @@ } /*}}}*/ /*FUNCTION Penta::GetPhi {{{*/ -void Penta::GetPhi(double* phi, double* epsilon, double viscosity){ +void Penta::GetPhi(IssmDouble* phi, IssmDouble* epsilon, IssmDouble viscosity){ /*Compute deformational heating from epsilon and viscosity */ - double epsilon_matrix[3][3]; - double epsilon_eff; - double epsilon_sqr[3][3]; + IssmDouble epsilon_matrix[3][3]; + IssmDouble epsilon_eff; + IssmDouble epsilon_sqr[3][3]; /* Build epsilon matrix */ epsilon_matrix[0][0]=*(epsilon+0); @@ -977,13 +977,13 @@ } /*}}}*/ /*FUNCTION Penta::GetStabilizationParameter {{{*/ -double Penta::GetStabilizationParameter(double u, double v, double w, double diameter, double kappa){ +IssmDouble Penta::GetStabilizationParameter(IssmDouble u, IssmDouble v, IssmDouble w, IssmDouble diameter, IssmDouble kappa){ /*Compute stabilization parameter*/ /*kappa=thermalconductivity/(rho_ice*hearcapacity) for thermal model*/ /*kappa=enthalpydiffusionparameter for enthalpy model*/ - double normu; - double tau_parameter; + IssmDouble normu; + IssmDouble tau_parameter; normu=pow(pow(u,2)+pow(v,2)+pow(w,2),0.5); if(normu*diameter/(3*2*kappa)<1){ @@ -995,7 +995,7 @@ } /*}}}*/ /*FUNCTION Penta::GetStrainRate3dPattyn{{{*/ -void Penta::GetStrainRate3dPattyn(double* epsilon,double* xyz_list, GaussPenta* gauss, Input* vx_input, Input* vy_input){ +void Penta::GetStrainRate3dPattyn(IssmDouble* epsilon,IssmDouble* xyz_list, GaussPenta* gauss, Input* vx_input, Input* vy_input){ /*Compute the 3d Blatter/PattynStrain Rate (5 components): * * epsilon=[exx eyy exy exz eyz] @@ -1007,8 +1007,8 @@ */ int i; - double epsilonvx[5]; - double epsilonvy[5]; + IssmDouble epsilonvx[5]; + IssmDouble epsilonvy[5]; /*Check that both inputs have been found*/ if (!vx_input || !vy_input){ @@ -1024,16 +1024,16 @@ } /*}}}*/ /*FUNCTION Penta::GetStrainRate3d{{{*/ -void Penta::GetStrainRate3d(double* epsilon,double* xyz_list, GaussPenta* gauss, Input* vx_input, Input* vy_input, Input* vz_input){ +void Penta::GetStrainRate3d(IssmDouble* epsilon,IssmDouble* xyz_list, GaussPenta* gauss, Input* vx_input, Input* vy_input, Input* vz_input){ /*Compute the 3d Strain Rate (6 components): * * epsilon=[exx eyy ezz exy exz eyz] */ int i; - double epsilonvx[6]; - double epsilonvy[6]; - double epsilonvz[6]; + IssmDouble epsilonvx[6]; + IssmDouble epsilonvy[6]; + IssmDouble epsilonvz[6]; /*Check that both inputs have been found*/ if (!vx_input || !vy_input || !vz_input){ @@ -1099,12 +1099,12 @@ } /*}}}*/ /*FUNCTION Penta::GetZcoord {{{*/ -double Penta::GetZcoord(GaussPenta* gauss){ +IssmDouble Penta::GetZcoord(GaussPenta* gauss){ int i; - double z; - double xyz_list[NUMVERTICES][3]; - double z_list[NUMVERTICES]; + IssmDouble z; + IssmDouble xyz_list[NUMVERTICES][3]; + IssmDouble z_list[NUMVERTICES]; GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES); for(i=0;iinputs->AddInput(new IntInput(name,(int)scalar)); } - else if ((code==7) || (code==3)){ //double - this->inputs->AddInput(new DoubleInput(name,(double)scalar)); + else if ((code==7) || (code==3)){ //IssmDouble + this->inputs->AddInput(new DoubleInput(name,(IssmDouble)scalar)); } else _error_("%s%i"," could not recognize nature of vector from code ",code); } /*}}}*/ -/*FUNCTION Penta::InputCreate(double* vector,int index,IoModel* iomodel,int M,int N,int vector_type,int vector_enum,int code){{{*/ -void Penta::InputCreate(double* vector, int index,IoModel* iomodel,int M,int N,int vector_type,int vector_enum,int code){//index into elements +/*FUNCTION Penta::InputCreate(IssmDouble* vector,int index,IoModel* iomodel,int M,int N,int vector_type,int vector_enum,int code){{{*/ +void Penta::InputCreate(IssmDouble* vector, int index,IoModel* iomodel,int M,int N,int vector_type,int vector_enum,int code){//index into elements /*Intermediaries*/ int i,j,t; int penta_vertex_ids[6]; int row; - double nodeinputs[6]; - double time; + IssmDouble nodeinputs[6]; + IssmDouble time; TransientInput* transientinput=NULL; int numberofvertices; int numberofelements; - double yts; + IssmDouble yts; /*Fetch parameters: */ iomodel->Constant(&numberofvertices,MeshNumberofverticesEnum); @@ -1221,7 +1221,7 @@ if(M==numberofvertices){ /*create input values: */ - for(i=0;i<6;i++)nodeinputs[i]=(double)vector[penta_vertex_ids[i]-1]; + for(i=0;i<6;i++)nodeinputs[i]=(IssmDouble)vector[penta_vertex_ids[i]-1]; /*process units: */ UnitConversion(&nodeinputs[0], 6 ,ExtToIuEnum, vector_enum); @@ -1236,14 +1236,14 @@ /*create input values: */ for(i=0;i<6;i++){ row=penta_vertex_ids[i]-1; - nodeinputs[i]=(double)vector[N*row+t]; + nodeinputs[i]=(IssmDouble)vector[N*row+t]; } /*process units: */ UnitConversion(&nodeinputs[0], 6 ,ExtToIuEnum, vector_enum); /*time? :*/ - time=(double)vector[(M-1)*N+t]*yts; + time=(IssmDouble)vector[(M-1)*N+t]*yts; if(t==0)transientinput=new TransientInput(vector_enum); transientinput->AddTimeInput(new PentaP1Input(vector_enum,nodeinputs),time); @@ -1264,8 +1264,8 @@ else if (code==6){ //integer this->inputs->AddInput(new IntInput(vector_enum,(int)vector[index])); } - else if (code==7){ //double - this->inputs->AddInput(new DoubleInput(vector_enum,(double)vector[index])); + else if (code==7){ //IssmDouble + this->inputs->AddInput(new DoubleInput(vector_enum,(IssmDouble)vector[index])); } else _error_("%s%i"," could not recognize nature of vector from code ",code); } @@ -1283,9 +1283,9 @@ void Penta::InputDepthAverageAtBase(int enum_type,int average_enum_type,int object_enum){ int step,i; - double xyz_list[NUMVERTICES][3]; - double Helem_list[NUMVERTICES]; - double zeros_list[NUMVERTICES]={0.0}; + IssmDouble xyz_list[NUMVERTICES][3]; + IssmDouble Helem_list[NUMVERTICES]; + IssmDouble zeros_list[NUMVERTICES]={0.0}; Penta* penta=NULL; Input* original_input=NULL; Input* element_integrated_input=NULL; @@ -1457,7 +1457,7 @@ } /*}}}*/ /*FUNCTION Penta::InputScale{{{*/ -void Penta::InputScale(int enum_type,double scale_factor){ +void Penta::InputScale(int enum_type,IssmDouble scale_factor){ Input* input=NULL; @@ -1470,7 +1470,7 @@ } /*}}}*/ /*FUNCTION Penta::InputToResult{{{*/ -void Penta::InputToResult(int enum_type,int step,double time){ +void Penta::InputToResult(int enum_type,int step,IssmDouble time){ int i; bool found = false; @@ -1504,8 +1504,8 @@ this->inputs->AddInput(new BoolInput(name,constant)); } /*}}}*/ -/*FUNCTION Penta::InputUpdateFromConstant(double value, int name);{{{*/ -void Penta::InputUpdateFromConstant(double constant, int name){ +/*FUNCTION Penta::InputUpdateFromConstant(IssmDouble value, int name);{{{*/ +void Penta::InputUpdateFromConstant(IssmDouble constant, int name){ /*Check that name is an element input*/ if (!IsInput(name)) return; @@ -1528,11 +1528,11 @@ /*Intermediaries*/ IssmInt i,j; int penta_vertex_ids[6]; - double nodeinputs[6]; - double cmmininputs[6]; - double cmmaxinputs[6]; + IssmDouble nodeinputs[6]; + IssmDouble cmmininputs[6]; + IssmDouble cmmaxinputs[6]; - double yts; + IssmDouble yts; bool control_analysis; int num_control_type; int num_cm_responses; @@ -1646,7 +1646,7 @@ } /*}}}*/ /*FUNCTION Penta::InputUpdateFromSolution {{{*/ -void Penta::InputUpdateFromSolution(double* solution){ +void Penta::InputUpdateFromSolution(IssmDouble* solution){ int analysis_type; @@ -1715,20 +1715,20 @@ } /*}}}*/ /*FUNCTION Penta::InputUpdateFromSolutionPrognostic{{{*/ -void Penta::InputUpdateFromSolutionPrognostic(double* solution){ +void Penta::InputUpdateFromSolutionPrognostic(IssmDouble* solution){ const int numdof = NDOF1*NUMVERTICES; const int numdof2d = NDOF1*NUMVERTICES2D; int i,hydroadjustment; int* doflist = NULL; - double rho_ice,rho_water,minthickness; - double newthickness[numdof]; - double newbed[numdof]; - double newsurface[numdof]; - double oldbed[NUMVERTICES]; - double oldsurface[NUMVERTICES]; - double oldthickness[NUMVERTICES]; + IssmDouble rho_ice,rho_water,minthickness; + IssmDouble newthickness[numdof]; + IssmDouble newbed[numdof]; + IssmDouble newsurface[numdof]; + IssmDouble oldbed[NUMVERTICES]; + IssmDouble oldsurface[NUMVERTICES]; + IssmDouble oldthickness[NUMVERTICES]; Penta *penta = NULL; /*If not on bed, return*/ @@ -1798,11 +1798,11 @@ } /*}}}*/ /*FUNCTION Penta::InputUpdateFromSolutionOneDof{{{*/ -void Penta::InputUpdateFromSolutionOneDof(double* solution,int enum_type){ +void Penta::InputUpdateFromSolutionOneDof(IssmDouble* solution,int enum_type){ const int numdof = NDOF1*NUMVERTICES; - double values[numdof]; + IssmDouble values[numdof]; int* doflist=NULL; /*Get dof list: */ @@ -1822,12 +1822,12 @@ } /*}}}*/ /*FUNCTION Penta::InputUpdateFromSolutionOneDofCollpased{{{*/ -void Penta::InputUpdateFromSolutionOneDofCollapsed(double* solution,int enum_type){ +void Penta::InputUpdateFromSolutionOneDofCollapsed(IssmDouble* solution,int enum_type){ const int numdof = NDOF1*NUMVERTICES; const int numdof2d = NDOF1*NUMVERTICES2D; - double values[numdof]; + IssmDouble values[numdof]; int* doflist = NULL; Penta *penta = NULL; @@ -1861,8 +1861,8 @@ xDelete(doflist); } /*}}}*/ -/*FUNCTION Penta::InputUpdateFromVector(double* vector, int name, int type);{{{*/ -void Penta::InputUpdateFromVector(double* vector, int name, int type){ +/*FUNCTION Penta::InputUpdateFromVector(IssmDouble* vector, int name, int type);{{{*/ +void Penta::InputUpdateFromVector(IssmDouble* vector, int name, int type){ /*Check that name is an element input*/ if (!IsInput(name)) return; @@ -1874,7 +1874,7 @@ case VertexEnum: /*New PentaVertexInpu*/ - double values[6]; + IssmDouble values[6]; /*Get values on the 6 vertices*/ for (int i=0;i<6;i++){ @@ -1977,7 +1977,7 @@ } /*}}}*/ /*FUNCTION Penta::IsNodeOnShelfFromFlags {{{*/ -bool Penta::IsNodeOnShelfFromFlags(double* flags){ +bool Penta::IsNodeOnShelfFromFlags(IssmDouble* flags){ int i; bool shelf=false; @@ -2008,14 +2008,14 @@ } /*}}}*/ /*FUNCTION Penta::ListResultsInfo{{{*/ -void Penta::ListResultsInfo(int** in_resultsenums,int** in_resultssizes,double** in_resultstimes,int** in_resultssteps,int* in_num_results){ +void Penta::ListResultsInfo(int** in_resultsenums,int** in_resultssizes,IssmDouble** in_resultstimes,int** in_resultssteps,int* in_num_results){ /*Intermediaries*/ int i; int numberofresults = 0; int *resultsenums = NULL; int *resultssizes = NULL; - double *resultstimes = NULL; + IssmDouble *resultstimes = NULL; int *resultssteps = NULL; /*Checks*/ @@ -2032,7 +2032,7 @@ /*Allocate output*/ resultsenums=xNew(numberofresults); resultssizes=xNew(numberofresults); - resultstimes=xNew(numberofresults); + resultstimes=xNew(numberofresults); resultssteps=xNew(numberofresults); /*populate enums*/ @@ -2059,14 +2059,14 @@ }/*}}}*/ /*FUNCTION Penta::MigrateGroundingLine{{{*/ -void Penta::MigrateGroundingLine(double* old_floating_ice,double* sheet_ungrounding){ +void Penta::MigrateGroundingLine(IssmDouble* old_floating_ice,IssmDouble* sheet_ungrounding){ int i,migration_style,unground; bool elementonshelf = false; - double bed_hydro,yts,gl_melting_rate; - double rho_water,rho_ice,density; - double melting[NUMVERTICES]; - double h[NUMVERTICES],s[NUMVERTICES],b[NUMVERTICES],ba[NUMVERTICES]; + IssmDouble bed_hydro,yts,gl_melting_rate; + IssmDouble rho_water,rho_ice,density; + IssmDouble melting[NUMVERTICES]; + IssmDouble h[NUMVERTICES],s[NUMVERTICES],b[NUMVERTICES],ba[NUMVERTICES]; if(!IsOnBed()) return; @@ -2142,13 +2142,13 @@ } /*}}}*/ /*FUNCTION Penta::MinEdgeLength{{{*/ -double Penta::MinEdgeLength(double xyz_list[6][3]){ +IssmDouble Penta::MinEdgeLength(IssmDouble xyz_list[6][3]){ /*Return the minimum lenght of the nine egdes of the penta*/ int i,node0,node1; int edges[9][2]={{0,1},{0,2},{1,2},{3,4},{3,5},{4,5},{0,3},{1,4},{2,5}}; //list of the nine edges - double length; - double minlength=-1; + IssmDouble length; + IssmDouble minlength=-1; for(i=0;i<9;i++){ /*Find the two nodes for this edge*/ @@ -2170,11 +2170,11 @@ } /*}}}*/ /*FUNCTION Penta::NodalValue {{{*/ -int Penta::NodalValue(double* pvalue, int index, int natureofdataenum,bool process_units){ +int Penta::NodalValue(IssmDouble* pvalue, int index, int natureofdataenum,bool process_units){ int i; int found=0; - double value; + IssmDouble value; Input* data=NULL; GaussPenta* gauss=NULL; @@ -2252,60 +2252,60 @@ } /*}}}*/ /*FUNCTION Penta::PositiveDegreeDay{{{*/ -void Penta::PositiveDegreeDay(double* pdds,double* pds,double signorm){ +void Penta::PositiveDegreeDay(IssmDouble* pdds,IssmDouble* pds,IssmDouble signorm){ int i,iqj,imonth; - double agd[NUMVERTICES]; // surface and basal - double saccu[NUMVERTICES] = {0}; // yearly surface accumulation - double smelt[NUMVERTICES] = {0}; // yearly melt - double precrunoff[NUMVERTICES]; // yearly runoff - double prect; // total precipitation during 1 year taking into account des. ef. - double water; //water=rain + snowmelt - double runoff; //meltwater only, does not include rain - double sconv; //rhow_rain/rhoi / 12 months + IssmDouble agd[NUMVERTICES]; // surface and basal + IssmDouble saccu[NUMVERTICES] = {0}; // yearly surface accumulation + IssmDouble smelt[NUMVERTICES] = {0}; // yearly melt + IssmDouble precrunoff[NUMVERTICES]; // yearly runoff + IssmDouble prect; // total precipitation during 1 year taking into account des. ef. + IssmDouble water; //water=rain + snowmelt + IssmDouble runoff; //meltwater only, does not include rain + IssmDouble sconv; //rhow_rain/rhoi / 12 months - double rho_water,rho_ice,density; - double lapser=6.5/1000, sealev=0; // lapse rate. degrees per meter. - double desfac = 0.5; //desert elevation factor - double s0p[NUMVERTICES]={0}; //should be set to elevation from precip source - double s0t[NUMVERTICES]={0}; //should be set to elevation from temperature source - double st; // elevation between altitude of the temp record and current altitude - double sp; // elevation between altitude of the prec record and current altitude + IssmDouble rho_water,rho_ice,density; + IssmDouble lapser=6.5/1000, sealev=0; // lapse rate. degrees per meter. + IssmDouble desfac = 0.5; //desert elevation factor + IssmDouble s0p[NUMVERTICES]={0}; //should be set to elevation from precip source + IssmDouble s0t[NUMVERTICES]={0}; //should be set to elevation from temperature source + IssmDouble st; // elevation between altitude of the temp record and current altitude + IssmDouble sp; // elevation between altitude of the prec record and current altitude // PDD and PD constants and variables - double siglim; // sigma limit for the integration which is equal to 2.5 sigmanorm - double signormc = signorm - 0.5; // sigma of the temperature distribution for cloudy day - double siglimc, siglim0, siglim0c; - double PDup, pddsig, PDCUT = 2.0; // PDcut: rain/snow cutoff temperature (C) - double DT = 0.02; - double pddt, pd; // pd: snow/precip fraction, precipitation falling as snow + IssmDouble siglim; // sigma limit for the integration which is equal to 2.5 sigmanorm + IssmDouble signormc = signorm - 0.5; // sigma of the temperature distribution for cloudy day + IssmDouble siglimc, siglim0, siglim0c; + IssmDouble PDup, pddsig, PDCUT = 2.0; // PDcut: rain/snow cutoff temperature (C) + IssmDouble DT = 0.02; + IssmDouble pddt, pd; // pd: snow/precip fraction, precipitation falling as snow - double q, qmpt; // q is desert/elev. fact, hnpfac is huybrect fact, and pd is normal dist. - double qm[NUMVERTICES] = {0}; // snow part of the precipitation - double qmt[NUMVERTICES] = {0}; // precipitation without desertification effect adjustment - double qmp[NUMVERTICES] = {0}; // desertification taken into account - double pdd[NUMVERTICES] = {0}; - double frzndd[NUMVERTICES] = {0}; + IssmDouble q, qmpt; // q is desert/elev. fact, hnpfac is huybrect fact, and pd is normal dist. + IssmDouble qm[NUMVERTICES] = {0}; // snow part of the precipitation + IssmDouble qmt[NUMVERTICES] = {0}; // precipitation without desertification effect adjustment + IssmDouble qmp[NUMVERTICES] = {0}; // desertification taken into account + IssmDouble pdd[NUMVERTICES] = {0}; + IssmDouble frzndd[NUMVERTICES] = {0}; - double tstar; // monthly mean surface temp - double Tsum[NUMVERTICES]= {0}; // average summer (JJA) temperature - double Tsurf[NUMVERTICES] = {0}; // average annual temperature + IssmDouble tstar; // monthly mean surface temp + IssmDouble Tsum[NUMVERTICES]= {0}; // average summer (JJA) temperature + IssmDouble Tsurf[NUMVERTICES] = {0}; // average annual temperature - double h[NUMVERTICES],s[NUMVERTICES],ttmp[NUMVERTICES],prectmp[NUMVERTICES]; // ,b[NUMVERTICES] - double t[NUMVERTICES][12],prec[NUMVERTICES][12]; - double deltm=1/12; + IssmDouble h[NUMVERTICES],s[NUMVERTICES],ttmp[NUMVERTICES],prectmp[NUMVERTICES]; // ,b[NUMVERTICES] + IssmDouble t[NUMVERTICES][12],prec[NUMVERTICES][12]; + IssmDouble deltm=1/12; int ismon[12]={12,1,2,3,4,5,6,7,8,9,10,11}; - double snwm; // snow that could have been melted in a year. - double snwmf; // ablation factor for snow per positive degree day. - double smf; // ablation factor for ice per pdd (Braithwaite 1995 from tarasov 2002). + IssmDouble snwm; // snow that could have been melted in a year. + IssmDouble snwmf; // ablation factor for snow per positive degree day. + IssmDouble smf; // ablation factor for ice per pdd (Braithwaite 1995 from tarasov 2002). - double dfrz=1.5, CovrLm=2009./3.35e+5, dCovrLm=dfrz*CovrLm; //m*J kg^-1 C^-1 /(J kg^-1)=m/C yr - double supice,supcap,diffndd; - double fsupT=0.5, fsupndd=0.5; // Tsurf mode factors for supice - double pddtj[NUMVERTICES], hmx2; + IssmDouble dfrz=1.5, CovrLm=2009./3.35e+5, dCovrLm=dfrz*CovrLm; //m*J kg^-1 C^-1 /(J kg^-1)=m/C yr + IssmDouble supice,supcap,diffndd; + IssmDouble fsupT=0.5, fsupndd=0.5; // Tsurf mode factors for supice + IssmDouble pddtj[NUMVERTICES], hmx2; /*Recover info at the vertices: */ GetInputListOnVertices(&h[0],ThicknessEnum); @@ -2480,9 +2480,9 @@ void Penta::PotentialSheetUngrounding(Vector* potential_sheet_ungrounding){ int i; - double h[NUMVERTICES],ba[NUMVERTICES]; - double bed_hydro; - double rho_water,rho_ice,density; + IssmDouble h[NUMVERTICES],ba[NUMVERTICES]; + IssmDouble bed_hydro; + IssmDouble rho_water,rho_ice,density; bool elementonshelf = false; /*material parameters: */ @@ -2517,15 +2517,15 @@ } /*}}}*/ /*FUNCTION Penta::ReduceMatrixStokes {{{*/ -void Penta::ReduceMatrixStokes(double* Ke_reduced, double* Ke_temp){ +void Penta::ReduceMatrixStokes(IssmDouble* Ke_reduced, IssmDouble* Ke_temp){ int i,j; - double Kii[24][24]; - double Kib[24][3]; - double Kbb[3][3]; - double Kbi[3][24]; - double Kbbinv[3][3]; - double Kright[24][24]; + IssmDouble Kii[24][24]; + IssmDouble Kib[24][3]; + IssmDouble Kbb[3][3]; + IssmDouble Kbi[3][24]; + IssmDouble Kbbinv[3][3]; + IssmDouble Kright[24][24]; /*Create the four matrices used for reduction */ for(i=0;i<24;i++){ @@ -2559,15 +2559,15 @@ } /*}}}*/ /*FUNCTION Penta::ReduceVectorStokes {{{*/ -void Penta::ReduceVectorStokes(double* Pe_reduced, double* Ke_temp, double* Pe_temp){ +void Penta::ReduceVectorStokes(IssmDouble* Pe_reduced, IssmDouble* Ke_temp, IssmDouble* Pe_temp){ int i,j; - double Pi[24]; - double Pb[3]; - double Kbb[3][3]; - double Kib[24][3]; - double Kbbinv[3][3]; - double Pright[24]; + IssmDouble Pi[24]; + IssmDouble Pb[3]; + IssmDouble Kbb[3][3]; + IssmDouble Kib[24][3]; + IssmDouble Kbbinv[3][3]; + IssmDouble Pright[24]; /*Create the four matrices used for reduction */ for(i=0;i<24;i++) Pi[i]=*(Pe_temp+i); @@ -2594,7 +2594,7 @@ } /*}}}*/ /*FUNCTION Penta::RequestedOutput{{{*/ -void Penta::RequestedOutput(int output_enum,int step,double time){ +void Penta::RequestedOutput(int output_enum,int step,IssmDouble time){ if(IsInput(output_enum)){ /*just transfer this input to results, and we are done: */ @@ -2648,9 +2648,9 @@ void Penta::ResetCoordinateSystem(void){ int approximation; - double slopex[NUMVERTICES]; - double slopey[NUMVERTICES]; - double xz_plane[6]; + IssmDouble slopex[NUMVERTICES]; + IssmDouble slopey[NUMVERTICES]; + IssmDouble xz_plane[6]; /*For Stokes only: we want the CS to be tangential to the bedrock*/ inputs->GetInputValue(&approximation,ApproximationEnum); @@ -2738,10 +2738,10 @@ } /*}}}*/ /*FUNCTION Penta::SurfaceArea {{{*/ -double Penta::SurfaceArea(void){ +IssmDouble Penta::SurfaceArea(void){ int approximation; - double S; + IssmDouble S; Tria* tria=NULL; /*retrieve inputs :*/ @@ -2775,12 +2775,12 @@ } /*}}}*/ /*FUNCTION Penta::SurfaceNormal {{{*/ -void Penta::SurfaceNormal(double* surface_normal, double xyz_list[3][3]){ +void Penta::SurfaceNormal(IssmDouble* surface_normal, IssmDouble xyz_list[3][3]){ int i; - double v13[3],v23[3]; - double normal[3]; - double normal_norm; + IssmDouble v13[3],v23[3]; + IssmDouble normal[3]; + IssmDouble normal_norm; for (i=0;i<3;i++){ v13[i]=xyz_list[0][i]-xyz_list[2][i]; @@ -2799,13 +2799,13 @@ } /*}}}*/ /*FUNCTION Penta::TimeAdapt{{{*/ -double Penta::TimeAdapt(void){ +IssmDouble Penta::TimeAdapt(void){ int i; - double C,dx,dy,dz,dt; - double maxabsvx,maxabsvy,maxabsvz; - double maxx,minx,maxy,miny,maxz,minz; - double xyz_list[NUMVERTICES][3]; + IssmDouble C,dx,dy,dz,dt; + IssmDouble maxabsvx,maxabsvy,maxabsvz; + IssmDouble maxx,minx,maxy,miny,maxz,minz; + IssmDouble xyz_list[NUMVERTICES][3]; /*get CFL coefficient:*/ this->parameters->FindParam(&C,TimesteppingCflCoefficientEnum); @@ -2850,12 +2850,12 @@ int penta_type; int penta_node_ids[6]; int penta_vertex_ids[6]; - double nodeinputs[6]; - double yts; + IssmDouble nodeinputs[6]; + IssmDouble yts; int stabilization; bool dakota_analysis; bool isstokes; - double beta,heatcapacity,referencetemperature,meltingpoint,latentheat; + IssmDouble beta,heatcapacity,referencetemperature,meltingpoint,latentheat; /*Fetch parameters: */ iomodel->Constant(&yts,ConstantsYtsEnum); @@ -2997,7 +2997,7 @@ } /*}}}*/ /*FUNCTION Penta::UpdatePotentialSheetUngrounding{{{*/ -int Penta::UpdatePotentialSheetUngrounding(double* vertices_potentially_ungrounding,Vector* vec_nodes_on_iceshelf,double* nodes_on_iceshelf){ +int Penta::UpdatePotentialSheetUngrounding(IssmDouble* vertices_potentially_ungrounding,Vector* vec_nodes_on_iceshelf,IssmDouble* nodes_on_iceshelf){ int i; int nflipped=0; @@ -3024,12 +3024,12 @@ /*Intermediaries*/ int iv; - double phi; - double viscosity; - double xyz_list[NUMVERTICES][3]; - double epsilon[6]; - double viscousheating[NUMVERTICES]={0,0,0,0,0,0}; - double thickness; + IssmDouble phi; + IssmDouble viscosity; + IssmDouble xyz_list[NUMVERTICES][3]; + IssmDouble epsilon[6]; + IssmDouble viscousheating[NUMVERTICES]={0,0,0,0,0,0}; + IssmDouble thickness; GaussPenta *gauss=NULL; /*Initialize Element vector*/ @@ -3065,18 +3065,18 @@ } /*}}}*/ /*FUNCTION Penta::SmearFunction {{{*/ -void Penta::SmearFunction(Vector* smearedvector,double (*WeightFunction)(double distance,double radius),double radius){ +void Penta::SmearFunction(Vector* smearedvector,IssmDouble (*WeightFunction)(IssmDouble distance,IssmDouble radius),IssmDouble radius){ _error_("not implemented yet"); } /*}}}*/ #ifdef _HAVE_RESPONSES_ /*FUNCTION Penta::IceVolume {{{*/ -double Penta::IceVolume(void){ +IssmDouble Penta::IceVolume(void){ /*The volume of a troncated prism is base * 1/3 sum(length of edges)*/ - double base,height; - double xyz_list[NUMVERTICES][3]; + IssmDouble base,height; + IssmDouble xyz_list[NUMVERTICES][3]; if(IsOnWater())return 0; @@ -3095,10 +3095,10 @@ } /*}}}*/ /*FUNCTION Penta::MinVel{{{*/ -void Penta::MinVel(double* pminvel, bool process_units){ +void Penta::MinVel(IssmDouble* pminvel, bool process_units){ /*Get minimum:*/ - double minvel=this->inputs->Min(VelEnum); + IssmDouble minvel=this->inputs->Min(VelEnum); /*process units if requested: */ if(process_units) minvel=UnitConversion(minvel,IuToExtEnum,VelEnum); @@ -3108,10 +3108,10 @@ } /*}}}*/ /*FUNCTION Penta::MinVx{{{*/ -void Penta::MinVx(double* pminvx, bool process_units){ +void Penta::MinVx(IssmDouble* pminvx, bool process_units){ /*Get minimum:*/ - double minvx=this->inputs->Min(VxEnum); + IssmDouble minvx=this->inputs->Min(VxEnum); /*process units if requested: */ if(process_units) minvx=UnitConversion(minvx,IuToExtEnum,VxEnum); @@ -3121,10 +3121,10 @@ } /*}}}*/ /*FUNCTION Penta::MinVy{{{*/ -void Penta::MinVy(double* pminvy, bool process_units){ +void Penta::MinVy(IssmDouble* pminvy, bool process_units){ /*Get minimum:*/ - double minvy=this->inputs->Min(VyEnum); + IssmDouble minvy=this->inputs->Min(VyEnum); /*process units if requested: */ if(process_units) minvy=UnitConversion(minvy,IuToExtEnum,VyEnum); @@ -3134,10 +3134,10 @@ } /*}}}*/ /*FUNCTION Penta::MinVz{{{*/ -void Penta::MinVz(double* pminvz, bool process_units){ +void Penta::MinVz(IssmDouble* pminvz, bool process_units){ /*Get minimum:*/ - double minvz=this->inputs->Min(VzEnum); + IssmDouble minvz=this->inputs->Min(VzEnum); /*process units if requested: */ if(process_units) minvz=UnitConversion(minvz,IuToExtEnum,VzEnum); @@ -3147,9 +3147,9 @@ } /*}}}*/ /*FUNCTION Penta::MassFlux {{{*/ -double Penta::MassFlux( double* segment,bool process_units){ +IssmDouble Penta::MassFlux( IssmDouble* segment,bool process_units){ - double mass_flux=0; + IssmDouble mass_flux=0; if(!IsOnBed()) return mass_flux; @@ -3171,10 +3171,10 @@ } /*}}}*/ /*FUNCTION Penta::MaxAbsVx{{{*/ -void Penta::MaxAbsVx(double* pmaxabsvx, bool process_units){ +void Penta::MaxAbsVx(IssmDouble* pmaxabsvx, bool process_units){ /*Get maximum:*/ - double maxabsvx=this->inputs->MaxAbs(VxEnum); + IssmDouble maxabsvx=this->inputs->MaxAbs(VxEnum); /*process units if requested: */ if(process_units) maxabsvx=UnitConversion(maxabsvx,IuToExtEnum,VxEnum); @@ -3184,10 +3184,10 @@ } /*}}}*/ /*FUNCTION Penta::MaxAbsVy{{{*/ -void Penta::MaxAbsVy(double* pmaxabsvy, bool process_units){ +void Penta::MaxAbsVy(IssmDouble* pmaxabsvy, bool process_units){ /*Get maximum:*/ - double maxabsvy=this->inputs->MaxAbs(VyEnum); + IssmDouble maxabsvy=this->inputs->MaxAbs(VyEnum); /*process units if requested: */ if(process_units) maxabsvy=UnitConversion(maxabsvy,IuToExtEnum,VyEnum); @@ -3197,10 +3197,10 @@ } /*}}}*/ /*FUNCTION Penta::MaxAbsVz{{{*/ -void Penta::MaxAbsVz(double* pmaxabsvz, bool process_units){ +void Penta::MaxAbsVz(IssmDouble* pmaxabsvz, bool process_units){ /*Get maximum:*/ - double maxabsvz=this->inputs->MaxAbs(VzEnum); + IssmDouble maxabsvz=this->inputs->MaxAbs(VzEnum); /*process units if requested: */ if(process_units) maxabsvz=UnitConversion(maxabsvz,IuToExtEnum,VyEnum); @@ -3210,10 +3210,10 @@ } /*}}}*/ /*FUNCTION Penta::MaxVel{{{*/ -void Penta::MaxVel(double* pmaxvel, bool process_units){ +void Penta::MaxVel(IssmDouble* pmaxvel, bool process_units){ /*Get maximum:*/ - double maxvel=this->inputs->Max(VelEnum); + IssmDouble maxvel=this->inputs->Max(VelEnum); /*process units if requested: */ if(process_units) maxvel=UnitConversion(maxvel,IuToExtEnum,VelEnum); @@ -3224,10 +3224,10 @@ } /*}}}*/ /*FUNCTION Penta::MaxVx{{{*/ -void Penta::MaxVx(double* pmaxvx, bool process_units){ +void Penta::MaxVx(IssmDouble* pmaxvx, bool process_units){ /*Get maximum:*/ - double maxvx=this->inputs->Max(VxEnum); + IssmDouble maxvx=this->inputs->Max(VxEnum); /*process units if requested: */ if(process_units) maxvx=UnitConversion(maxvx,IuToExtEnum,VxEnum); @@ -3237,10 +3237,10 @@ } /*}}}*/ /*FUNCTION Penta::MaxVy{{{*/ -void Penta::MaxVy(double* pmaxvy, bool process_units){ +void Penta::MaxVy(IssmDouble* pmaxvy, bool process_units){ /*Get maximum:*/ - double maxvy=this->inputs->Max(VyEnum); + IssmDouble maxvy=this->inputs->Max(VyEnum); /*process units if requested: */ if(process_units) maxvy=UnitConversion(maxvy,IuToExtEnum,VyEnum); @@ -3250,10 +3250,10 @@ } /*}}}*/ /*FUNCTION Penta::MaxVz{{{*/ -void Penta::MaxVz(double* pmaxvz, bool process_units){ +void Penta::MaxVz(IssmDouble* pmaxvz, bool process_units){ /*Get maximum:*/ - double maxvz=this->inputs->Max(VzEnum); + IssmDouble maxvz=this->inputs->Max(VzEnum); /*process units if requested: */ if(process_units) maxvz=UnitConversion(maxvz,IuToExtEnum,VzEnum); @@ -3263,7 +3263,7 @@ } /*}}}*/ /*FUNCTION Penta::ElementResponse{{{*/ -void Penta::ElementResponse(double* presponse,int response_enum,bool process_units){ +void Penta::ElementResponse(IssmDouble* presponse,int response_enum,bool process_units){ switch(response_enum){ case MaterialsRheologyBbarEnum: @@ -3272,7 +3272,7 @@ case VelEnum: /*Get input:*/ - double vel; + IssmDouble vel; Input* vel_input; vel_input=this->inputs->GetInput(VelEnum); _assert_(vel_input); @@ -3315,24 +3315,24 @@ /*Intermediaries */ int stabilization; int i,j,ig,found=0; - double Jdet,u,v,w,um,vm,wm; - double h,hx,hy,hz,vx,vy,vz,vel; - double gravity,rho_ice,rho_water; - double epsvel=2.220446049250313e-16; - double heatcapacity,thermalconductivity,dt; - double pressure,enthalpy; - double latentheat,kappa; - double tau_parameter,diameter; - double xyz_list[NUMVERTICES][3]; - double B_conduct[3][numdof]; - double B_advec[3][numdof]; - double Bprime_advec[3][numdof]; - double L[numdof]; - double dbasis[3][6]; - double D_scalar_conduct,D_scalar_advec; - double D_scalar_trans,D_scalar_stab; - double D[3][3]; - double K[3][3]={0.0}; + IssmDouble Jdet,u,v,w,um,vm,wm; + IssmDouble h,hx,hy,hz,vx,vy,vz,vel; + IssmDouble gravity,rho_ice,rho_water; + IssmDouble epsvel=2.220446049250313e-16; + IssmDouble heatcapacity,thermalconductivity,dt; + IssmDouble pressure,enthalpy; + IssmDouble latentheat,kappa; + IssmDouble tau_parameter,diameter; + IssmDouble xyz_list[NUMVERTICES][3]; + IssmDouble B_conduct[3][numdof]; + IssmDouble B_advec[3][numdof]; + IssmDouble Bprime_advec[3][numdof]; + IssmDouble L[numdof]; + IssmDouble dbasis[3][6]; + IssmDouble D_scalar_conduct,D_scalar_advec; + IssmDouble D_scalar_trans,D_scalar_stab; + IssmDouble D[3][3]; + IssmDouble K[3][3]={0.0}; Tria* tria=NULL; GaussPenta *gauss=NULL; @@ -3471,13 +3471,13 @@ /*Intermediaries */ int i,j,ig; - double mixed_layer_capacity,thermal_exchange_velocity; - double rho_ice,rho_water,heatcapacity; - double Jdet2d,dt; - double xyz_list[NUMVERTICES][3]; - double xyz_list_tria[NUMVERTICES2D][3]; - double basis[NUMVERTICES]; - double D_scalar; + IssmDouble mixed_layer_capacity,thermal_exchange_velocity; + IssmDouble rho_ice,rho_water,heatcapacity; + IssmDouble Jdet2d,dt; + IssmDouble xyz_list[NUMVERTICES][3]; + IssmDouble xyz_list_tria[NUMVERTICES2D][3]; + IssmDouble basis[NUMVERTICES]; + IssmDouble D_scalar; GaussPenta *gauss=NULL; /*Initialize Element matrix and return if necessary*/ @@ -3552,21 +3552,21 @@ /*Intermediaries */ int stabilization; int i,j,ig,found=0; - double Jdet,u,v,w,um,vm,wm,vel; - double h,hx,hy,hz,vx,vy,vz; - double gravity,rho_ice,rho_water,kappa; - double heatcapacity,thermalconductivity,dt; - double tau_parameter,diameter; - double xyz_list[NUMVERTICES][3]; - double B_conduct[3][numdof]; - double B_advec[3][numdof]; - double Bprime_advec[3][numdof]; - double L[numdof]; - double dbasis[3][6]; - double D_scalar_conduct,D_scalar_advec; - double D_scalar_trans,D_scalar_stab; - double D[3][3]; - double K[3][3]={0.0}; + IssmDouble Jdet,u,v,w,um,vm,wm,vel; + IssmDouble h,hx,hy,hz,vx,vy,vz; + IssmDouble gravity,rho_ice,rho_water,kappa; + IssmDouble heatcapacity,thermalconductivity,dt; + IssmDouble tau_parameter,diameter; + IssmDouble xyz_list[NUMVERTICES][3]; + IssmDouble B_conduct[3][numdof]; + IssmDouble B_advec[3][numdof]; + IssmDouble Bprime_advec[3][numdof]; + IssmDouble L[numdof]; + IssmDouble dbasis[3][6]; + IssmDouble D_scalar_conduct,D_scalar_advec; + IssmDouble D_scalar_trans,D_scalar_stab; + IssmDouble D[3][3]; + IssmDouble K[3][3]={0.0}; Tria* tria=NULL; GaussPenta *gauss=NULL; @@ -3704,13 +3704,13 @@ /*Intermediaries */ int i,j,ig; - double mixed_layer_capacity,thermal_exchange_velocity; - double rho_ice,rho_water,heatcapacity; - double Jdet2d,dt; - double xyz_list[NUMVERTICES][3]; - double xyz_list_tria[NUMVERTICES2D][3]; - double basis[NUMVERTICES]; - double D_scalar; + IssmDouble mixed_layer_capacity,thermal_exchange_velocity; + IssmDouble rho_ice,rho_water,heatcapacity; + IssmDouble Jdet2d,dt; + IssmDouble xyz_list[NUMVERTICES][3]; + IssmDouble xyz_list_tria[NUMVERTICES2D][3]; + IssmDouble basis[NUMVERTICES]; + IssmDouble D_scalar; GaussPenta *gauss=NULL; /*Initialize Element matrix and return if necessary*/ @@ -3775,19 +3775,19 @@ /*Intermediaries*/ int i,j,ig,found=0; int friction_type,stabilization; - double Jdet,phi,dt; - double rho_ice,heatcapacity; - double thermalconductivity,kappa; - double viscosity,pressure; - double enthalpy,enthalpypicard; - double tau_parameter,diameter; - double u,v,w; - double scalar_def,scalar_transient; - double temperature_list[NUMVERTICES]; - double xyz_list[NUMVERTICES][3]; - double L[numdof]; - double dbasis[3][6]; - double epsilon[6]; + IssmDouble Jdet,phi,dt; + IssmDouble rho_ice,heatcapacity; + IssmDouble thermalconductivity,kappa; + IssmDouble viscosity,pressure; + IssmDouble enthalpy,enthalpypicard; + IssmDouble tau_parameter,diameter; + IssmDouble u,v,w; + IssmDouble scalar_def,scalar_transient; + IssmDouble temperature_list[NUMVERTICES]; + IssmDouble xyz_list[NUMVERTICES][3]; + IssmDouble L[numdof]; + IssmDouble dbasis[3][6]; + IssmDouble epsilon[6]; GaussPenta *gauss=NULL; /*Initialize Element vector*/ @@ -3870,13 +3870,13 @@ /*Intermediaries */ int i,j,ig; - double Jdet2d; - double heatcapacity,h_pmp; - double mixed_layer_capacity,thermal_exchange_velocity; - double rho_ice,rho_water,pressure,dt,scalar_ocean; - double xyz_list[NUMVERTICES][3]; - double xyz_list_tria[NUMVERTICES2D][3]; - double basis[NUMVERTICES]; + IssmDouble Jdet2d; + IssmDouble heatcapacity,h_pmp; + IssmDouble mixed_layer_capacity,thermal_exchange_velocity; + IssmDouble rho_ice,rho_water,pressure,dt,scalar_ocean; + IssmDouble xyz_list[NUMVERTICES][3]; + IssmDouble xyz_list_tria[NUMVERTICES2D][3]; + IssmDouble basis[NUMVERTICES]; GaussPenta* gauss=NULL; /* Ice/ocean heat exchange flux on ice shelf base */ @@ -3928,14 +3928,14 @@ /*Intermediaries */ int i,j,ig; int analysis_type; - double xyz_list[NUMVERTICES][3]; - double xyz_list_tria[NUMVERTICES2D][3]={0.0}; - double Jdet2d,dt; - double rho_ice,heatcapacity,geothermalflux_value; - double basalfriction,alpha2,vx,vy; - double scalar,enthalpy,enthalpyup; - double pressure,pressureup; - double basis[NUMVERTICES]; + IssmDouble xyz_list[NUMVERTICES][3]; + IssmDouble xyz_list_tria[NUMVERTICES2D][3]={0.0}; + IssmDouble Jdet2d,dt; + IssmDouble rho_ice,heatcapacity,geothermalflux_value; + IssmDouble basalfriction,alpha2,vx,vy; + IssmDouble scalar,enthalpy,enthalpyup; + IssmDouble pressure,pressureup; + IssmDouble basis[NUMVERTICES]; Friction* friction=NULL; GaussPenta* gauss=NULL; GaussPenta* gaussup=NULL; @@ -4038,18 +4038,18 @@ /*Intermediaries*/ int i,j,ig,found=0; int friction_type,stabilization; - double Jdet,phi,dt; - double rho_ice,heatcapacity; - double thermalconductivity,kappa; - double viscosity,temperature; - double tau_parameter,diameter; - double u,v,w; - double scalar_def,scalar_transient; - double temperature_list[NUMVERTICES]; - double xyz_list[NUMVERTICES][3]; - double L[numdof]; - double dbasis[3][6]; - double epsilon[6]; + IssmDouble Jdet,phi,dt; + IssmDouble rho_ice,heatcapacity; + IssmDouble thermalconductivity,kappa; + IssmDouble viscosity,temperature; + IssmDouble tau_parameter,diameter; + IssmDouble u,v,w; + IssmDouble scalar_def,scalar_transient; + IssmDouble temperature_list[NUMVERTICES]; + IssmDouble xyz_list[NUMVERTICES][3]; + IssmDouble L[numdof]; + IssmDouble dbasis[3][6]; + IssmDouble epsilon[6]; GaussPenta *gauss=NULL; /*Initialize Element vector*/ @@ -4124,13 +4124,13 @@ /*Intermediaries */ int i,j,ig; - double Jdet2d; - double mixed_layer_capacity,thermal_exchange_velocity; - double rho_ice,rho_water,pressure,dt,scalar_ocean; - double heatcapacity,t_pmp; - double xyz_list[NUMVERTICES][3]; - double xyz_list_tria[NUMVERTICES2D][3]; - double basis[NUMVERTICES]; + IssmDouble Jdet2d; + IssmDouble mixed_layer_capacity,thermal_exchange_velocity; + IssmDouble rho_ice,rho_water,pressure,dt,scalar_ocean; + IssmDouble heatcapacity,t_pmp; + IssmDouble xyz_list[NUMVERTICES][3]; + IssmDouble xyz_list_tria[NUMVERTICES2D][3]; + IssmDouble basis[NUMVERTICES]; GaussPenta* gauss=NULL; /* Ice/ocean heat exchange flux on ice shelf base */ @@ -4182,13 +4182,13 @@ /*Intermediaries */ int i,j,ig; int analysis_type; - double xyz_list[NUMVERTICES][3]; - double xyz_list_tria[NUMVERTICES2D][3]={0.0}; - double Jdet2d,dt; - double rho_ice,heatcapacity,geothermalflux_value; - double basalfriction,alpha2,vx,vy; - double basis[NUMVERTICES]; - double scalar; + IssmDouble xyz_list[NUMVERTICES][3]; + IssmDouble xyz_list_tria[NUMVERTICES2D][3]={0.0}; + IssmDouble Jdet2d,dt; + IssmDouble rho_ice,heatcapacity,geothermalflux_value; + IssmDouble basalfriction,alpha2,vx,vy; + IssmDouble basis[NUMVERTICES]; + IssmDouble scalar; Friction* friction=NULL; GaussPenta* gauss=NULL; @@ -4247,8 +4247,8 @@ int i; int* doflist=NULL; - double values[numdof]; - double temp; + IssmDouble values[numdof]; + IssmDouble temp; GaussPenta *gauss=NULL; /*Get dof list: */ @@ -4278,8 +4278,8 @@ int i; int* doflist=NULL; - double values[numdof]; - double enthalpy; + IssmDouble values[numdof]; + IssmDouble enthalpy; GaussPenta *gauss=NULL; /*Get dof list: */ @@ -4303,18 +4303,18 @@ } /*}}}*/ /*FUNCTION Penta::InputUpdateFromSolutionThermal {{{*/ -void Penta::InputUpdateFromSolutionThermal(double* solution){ +void Penta::InputUpdateFromSolutionThermal(IssmDouble* solution){ const int numdof=NDOF1*NUMVERTICES; bool converged; int i,rheology_law; - double xyz_list[NUMVERTICES][3]; - double values[numdof]; - double B[numdof]; - double B_average,s_average; + IssmDouble xyz_list[NUMVERTICES][3]; + IssmDouble values[numdof]; + IssmDouble B[numdof]; + IssmDouble B_average,s_average; int* doflist=NULL; - //double pressure[numdof]; + //IssmDouble pressure[numdof]; /*Get dof list: */ GetDofList(&doflist,NoneApproximationEnum,GsetEnum); @@ -4374,19 +4374,19 @@ } /*}}}*/ /*FUNCTION Penta::InputUpdateFromSolutionEnthalpy {{{*/ -void Penta::InputUpdateFromSolutionEnthalpy(double* solution){ +void Penta::InputUpdateFromSolutionEnthalpy(IssmDouble* solution){ const int numdof=NDOF1*NUMVERTICES; bool converged=false; int i,rheology_law; - double xyz_list[NUMVERTICES][3]; - double values[numdof]; - double pressure[NUMVERTICES]; - double temperatures[numdof]; - double waterfraction[numdof]; - double B[numdof]; - double B_average,s_average; + IssmDouble xyz_list[NUMVERTICES][3]; + IssmDouble values[numdof]; + IssmDouble pressure[NUMVERTICES]; + IssmDouble temperatures[numdof]; + IssmDouble waterfraction[numdof]; + IssmDouble B[numdof]; + IssmDouble B_average,s_average; int* doflist=NULL; /*Get dof list: */ @@ -4475,7 +4475,7 @@ }/*}}}*/ /*FUNCTION Penta::ControlInputScaleGradient{{{*/ -void Penta::ControlInputScaleGradient(int enum_type,double scale){ +void Penta::ControlInputScaleGradient(int enum_type,IssmDouble scale){ Input* input=NULL; @@ -4491,10 +4491,10 @@ ((ControlInput*)input)->ScaleGradient(scale); }/*}}}*/ /*FUNCTION Penta::ControlInputSetGradient{{{*/ -void Penta::ControlInputSetGradient(double* gradient,int enum_type,int control_index){ +void Penta::ControlInputSetGradient(IssmDouble* gradient,int enum_type,int control_index){ int doflist1[NUMVERTICES]; - double grad_list[NUMVERTICES]; + IssmDouble grad_list[NUMVERTICES]; Input* grad_input=NULL; Input* input=NULL; @@ -4565,15 +4565,15 @@ /*Intermediaries */ int i,j,ig; bool incomplete_adjoint; - double xyz_list[NUMVERTICES][3]; - double Jdet; - double eps1dotdphii,eps1dotdphij; - double eps2dotdphii,eps2dotdphij; - double mu_prime; - double epsilon[5]; /* epsilon=[exx,eyy,exy,exz,eyz];*/ - double eps1[3],eps2[3]; - double phi[NUMVERTICES]; - double dphi[3][NUMVERTICES]; + IssmDouble xyz_list[NUMVERTICES][3]; + IssmDouble Jdet; + IssmDouble eps1dotdphii,eps1dotdphij; + IssmDouble eps2dotdphii,eps2dotdphij; + IssmDouble mu_prime; + IssmDouble epsilon[5]; /* epsilon=[exx,eyy,exy,exz,eyz];*/ + IssmDouble eps1[3],eps2[3]; + IssmDouble phi[NUMVERTICES]; + IssmDouble dphi[3][NUMVERTICES]; GaussPenta *gauss=NULL; /*Initialize Jacobian with regular Pattyn (first part of the Gateau derivative)*/ @@ -4633,16 +4633,16 @@ /*Intermediaries */ int i,j,ig; bool incomplete_adjoint; - double xyz_list[NUMVERTICES][3]; - double Jdet; - double eps1dotdphii,eps1dotdphij; - double eps2dotdphii,eps2dotdphij; - double eps3dotdphii,eps3dotdphij; - double mu_prime; - double epsilon[5]; /* epsilon=[exx,eyy,exy,exz,eyz];*/ - double eps1[3],eps2[3],eps3[3]; - double phi[NUMVERTICES]; - double dphi[3][NUMVERTICES]; + IssmDouble xyz_list[NUMVERTICES][3]; + IssmDouble Jdet; + IssmDouble eps1dotdphii,eps1dotdphij; + IssmDouble eps2dotdphii,eps2dotdphij; + IssmDouble eps3dotdphii,eps3dotdphij; + IssmDouble mu_prime; + IssmDouble epsilon[5]; /* epsilon=[exx,eyy,exy,exz,eyz];*/ + IssmDouble eps1[3],eps2[3],eps3[3]; + IssmDouble phi[NUMVERTICES]; + IssmDouble dphi[3][NUMVERTICES]; GaussPenta *gauss=NULL; /*Initialize Jacobian with regular Stokes (first part of the Gateau derivative)*/ @@ -4888,14 +4888,14 @@ int i,j,ig; int analysis_type; int doflist1[NUMVERTICES]; - double vx,vy,lambda,mu,alpha_complement,Jdet; - double bed,thickness,Neff,drag; - double xyz_list[NUMVERTICES][3]; - double xyz_list_tria[NUMVERTICES2D][3]={0.0}; - double dk[NDOF3]; - double grade_g[NUMVERTICES]={0.0}; - double grade_g_gaussian[NUMVERTICES]; - double basis[6]; + IssmDouble vx,vy,lambda,mu,alpha_complement,Jdet; + IssmDouble bed,thickness,Neff,drag; + IssmDouble xyz_list[NUMVERTICES][3]; + IssmDouble xyz_list_tria[NUMVERTICES2D][3]={0.0}; + IssmDouble dk[NDOF3]; + IssmDouble grade_g[NUMVERTICES]={0.0}; + IssmDouble grade_g_gaussian[NUMVERTICES]; + IssmDouble basis[6]; Friction* friction=NULL; GaussPenta *gauss=NULL; @@ -4959,16 +4959,16 @@ int i,j,ig; int analysis_type; int doflist1[NUMVERTICES]; - double bed,thickness,Neff; - double lambda,mu,xi,Jdet,vx,vy,vz; - double alpha_complement,drag; - double surface_normal[3],bed_normal[3]; - double xyz_list[NUMVERTICES][3]; - double xyz_list_tria[NUMVERTICES2D][3]={0.0}; - double dk[NDOF3]; - double basis[6]; - double grade_g[NUMVERTICES]={0.0}; - double grade_g_gaussian[NUMVERTICES]; + IssmDouble bed,thickness,Neff; + IssmDouble lambda,mu,xi,Jdet,vx,vy,vz; + IssmDouble alpha_complement,drag; + IssmDouble surface_normal[3],bed_normal[3]; + IssmDouble xyz_list[NUMVERTICES][3]; + IssmDouble xyz_list_tria[NUMVERTICES2D][3]={0.0}; + IssmDouble dk[NDOF3]; + IssmDouble basis[6]; + IssmDouble grade_g[NUMVERTICES]={0.0}; + IssmDouble grade_g_gaussian[NUMVERTICES]; Friction* friction=NULL; GaussPenta* gauss=NULL; @@ -5098,7 +5098,7 @@ this->matice->inputs->DeleteInput(MaterialsRheologyBbarEnum); } /*}}}*/ /*FUNCTION Penta::InputControlUpdate{{{*/ -void Penta::InputControlUpdate(double scalar,bool save_parameter){ +void Penta::InputControlUpdate(IssmDouble scalar,bool save_parameter){ /*Intermediary*/ int num_controls; @@ -5136,16 +5136,16 @@ } /*}}}*/ /*FUNCTION Penta::InputUpdateFromSolutionAdjointStokes {{{*/ -void Penta::InputUpdateFromSolutionAdjointStokes(double* solution){ +void Penta::InputUpdateFromSolutionAdjointStokes(IssmDouble* solution){ const int numdof=NDOF4*NUMVERTICES; int i; - double values[numdof]; - double lambdax[NUMVERTICES]; - double lambday[NUMVERTICES]; - double lambdaz[NUMVERTICES]; - double lambdap[NUMVERTICES]; + IssmDouble values[numdof]; + IssmDouble lambdax[NUMVERTICES]; + IssmDouble lambday[NUMVERTICES]; + IssmDouble lambdaz[NUMVERTICES]; + IssmDouble lambdap[NUMVERTICES]; int* doflist=NULL; /*Get dof list: */ @@ -5179,14 +5179,14 @@ } /*}}}*/ /*FUNCTION Penta::InputUpdateFromSolutionAdjointHoriz {{{*/ -void Penta::InputUpdateFromSolutionAdjointHoriz(double* solution){ +void Penta::InputUpdateFromSolutionAdjointHoriz(IssmDouble* solution){ const int numdof=NDOF2*NUMVERTICES; int i; - double values[numdof]; - double lambdax[NUMVERTICES]; - double lambday[NUMVERTICES]; + IssmDouble values[numdof]; + IssmDouble lambdax[NUMVERTICES]; + IssmDouble lambday[NUMVERTICES]; int* doflist=NULL; /*Get dof list: */ @@ -5214,10 +5214,10 @@ } /*}}}*/ /*FUNCTION Penta::SurfaceAverageVelMisfit {{{*/ -double Penta::SurfaceAverageVelMisfit(bool process_units,int weight_index){ +IssmDouble Penta::SurfaceAverageVelMisfit(bool process_units,int weight_index){ int approximation; - double J; + IssmDouble J; Tria* tria=NULL; /*retrieve inputs :*/ @@ -5251,10 +5251,10 @@ } /*}}}*/ /*FUNCTION Penta::SurfaceAbsVelMisfit {{{*/ -double Penta::SurfaceAbsVelMisfit(bool process_units,int weight_index){ +IssmDouble Penta::SurfaceAbsVelMisfit(bool process_units,int weight_index){ int approximation; - double J; + IssmDouble J; Tria* tria=NULL; /*retrieve inputs :*/ @@ -5288,10 +5288,10 @@ } /*}}}*/ /*FUNCTION Penta::SurfaceLogVelMisfit {{{*/ -double Penta::SurfaceLogVelMisfit(bool process_units,int weight_index){ +IssmDouble Penta::SurfaceLogVelMisfit(bool process_units,int weight_index){ int approximation; - double J; + IssmDouble J; Tria* tria=NULL; /*retrieve inputs :*/ @@ -5325,9 +5325,9 @@ } /*}}}*/ /*FUNCTION Penta::SurfaceLogVxVyMisfit {{{*/ -double Penta::SurfaceLogVxVyMisfit(bool process_units,int weight_index){ +IssmDouble Penta::SurfaceLogVxVyMisfit(bool process_units,int weight_index){ - double J; + IssmDouble J; Tria* tria=NULL; /*inputs: */ @@ -5364,10 +5364,10 @@ } /*}}}*/ /*FUNCTION Penta::SurfaceRelVelMisfit {{{*/ -double Penta::SurfaceRelVelMisfit(bool process_units,int weight_index){ +IssmDouble Penta::SurfaceRelVelMisfit(bool process_units,int weight_index){ int approximation; - double J; + IssmDouble J; Tria* tria=NULL; /*retrieve inputs :*/ @@ -5401,16 +5401,16 @@ } /*}}}*/ /*FUNCTION Penta::ThicknessAbsGradient{{{*/ -double Penta::ThicknessAbsGradient(bool process_units,int weight_index){ +IssmDouble Penta::ThicknessAbsGradient(bool process_units,int weight_index){ _error_("Not implemented yet"); } /*}}}*/ /*FUNCTION Penta::ThicknessAbsMisfit {{{*/ -double Penta::ThicknessAbsMisfit(bool process_units,int weight_index){ +IssmDouble Penta::ThicknessAbsMisfit(bool process_units,int weight_index){ int approximation; - double J; + IssmDouble J; Tria* tria=NULL; /*retrieve inputs :*/ @@ -5427,9 +5427,9 @@ } /*}}}*/ /*FUNCTION Penta::DragCoefficientAbsGradient{{{*/ -double Penta::DragCoefficientAbsGradient(bool process_units,int weight_index){ +IssmDouble Penta::DragCoefficientAbsGradient(bool process_units,int weight_index){ - double J; + IssmDouble J; Tria* tria=NULL; /*If on water, on shelf or not on bed, skip: */ @@ -5442,9 +5442,9 @@ } /*}}}*/ /*FUNCTION Penta::RheologyBbarAbsGradient{{{*/ -double Penta::RheologyBbarAbsGradient(bool process_units,int weight_index){ +IssmDouble Penta::RheologyBbarAbsGradient(bool process_units,int weight_index){ - double J; + IssmDouble J; Tria* tria=NULL; /*If on water, on shelf or not on bed, skip: */ @@ -5480,9 +5480,9 @@ } /*}}}*/ /*FUNCTION Penta::SetControlInputsFromVector{{{*/ -void Penta::SetControlInputsFromVector(double* vector,int control_enum,int control_index){ +void Penta::SetControlInputsFromVector(IssmDouble* vector,int control_enum,int control_index){ - double values[NUMVERTICES]; + IssmDouble values[NUMVERTICES]; int doflist1[NUMVERTICES]; Input *input = NULL; Input *new_input = NULL; @@ -5517,8 +5517,8 @@ #endif #ifdef _HAVE_DAKOTA_ -/*FUNCTION Penta::InputUpdateFromVectorDakota(double* vector, int name, int type);{{{*/ -void Penta::InputUpdateFromVectorDakota(double* vector, int name, int type){ +/*FUNCTION Penta::InputUpdateFromVectorDakota(IssmDouble* vector, int name, int type);{{{*/ +void Penta::InputUpdateFromVectorDakota(IssmDouble* vector, int name, int type){ int i,j; @@ -5530,7 +5530,7 @@ case VertexEnum: /*New PentaP1Input*/ - double values[6]; + IssmDouble values[6]; /*Get values on the 6 vertices*/ for (i=0;i<6;i++){ @@ -5541,11 +5541,11 @@ switch(name){ case ThicknessEnum: /*Update thickness + surface: assume bed is constant. On ice shelves, takes hydrostatic equilibrium {{{*/ - double thickness[6]; - double thickness_init[6]; - double hydrostatic_ratio[6]; - double surface[6]; - double bed[6]; + IssmDouble thickness[6]; + IssmDouble thickness_init[6]; + IssmDouble hydrostatic_ratio[6]; + IssmDouble surface[6]; + IssmDouble bed[6]; /*retrieve inputs: */ GetInputListOnVertices(&thickness_init[0],ThicknessEnum); @@ -5559,7 +5559,7 @@ /*build new bed and surface: */ if (this->IsFloating()){ /*hydrostatic equilibrium: */ - double rho_ice,rho_water,di; + IssmDouble rho_ice,rho_water,di; rho_ice=this->matpar->GetRhoIce(); rho_water=this->matpar->GetRhoWater(); @@ -5629,15 +5629,15 @@ _error_(" not supported yet!"); } /*}}}*/ -/*FUNCTION Penta::InputUpdateFromMatrixDakota(double* matrix, int nrows, int ncols, int name, int type);{{{*/ -void Penta::InputUpdateFromMatrixDakota(double* matrix, int nrows, int ncols, int name, int type){ +/*FUNCTION Penta::InputUpdateFromMatrixDakota(IssmDouble* matrix, int nrows, int ncols, int name, int type);{{{*/ +void Penta::InputUpdateFromMatrixDakota(IssmDouble* matrix, int nrows, int ncols, int name, int type){ int i,j,t; TransientInput* transientinput=NULL; - double values[6]; - double time; + IssmDouble values[6]; + IssmDouble time; int row; - double yts; + IssmDouble yts; /*Check that name is an element input*/ if (!IsInput(name)) return; @@ -5655,11 +5655,11 @@ /*create input values: */ for(i=0;i<6;i++){ row=this->nodes[i]->GetSidList(); - values[i]=(double)matrix[ncols*row+t]; + values[i]=(IssmDouble)matrix[ncols*row+t]; } /*time? :*/ - time=(double)matrix[(nrows-1)*ncols+t]*yts; + time=(IssmDouble)matrix[(nrows-1)*ncols+t]*yts; if(t==0) transientinput=new TransientInput(name); transientinput->AddTimeInput(new PentaP1Input(name,values),time); @@ -5741,16 +5741,16 @@ /*Intermediaries */ int i,j,ig; - double Jdet; - double viscosity,oldviscosity,newviscosity,viscosity_overshoot; //viscosity - double epsilon[5],oldepsilon[5]; /* epsilon=[exx,eyy,exy,exz,eyz];*/ - double xyz_list[NUMVERTICES][3]; - double B[3][numdofp]; - double Bprime[3][numdofm]; - double D[3][3]={0.0}; // material matrix, simple scalar matrix. - double D_scalar; - double Ke_gg[numdofp][numdofm]={0.0}; //local element stiffness matrix - double Ke_gg_gaussian[numdofp][numdofm]; //stiffness matrix evaluated at the gaussian point. + IssmDouble Jdet; + IssmDouble viscosity,oldviscosity,newviscosity,viscosity_overshoot; //viscosity + IssmDouble epsilon[5],oldepsilon[5]; /* epsilon=[exx,eyy,exy,exz,eyz];*/ + IssmDouble xyz_list[NUMVERTICES][3]; + IssmDouble B[3][numdofp]; + IssmDouble Bprime[3][numdofm]; + IssmDouble D[3][3]={0.0}; // material matrix, simple scalar matrix. + IssmDouble D_scalar; + IssmDouble Ke_gg[numdofp][numdofm]={0.0}; //local element stiffness matrix + IssmDouble Ke_gg_gaussian[numdofp][numdofm]; //stiffness matrix evaluated at the gaussian point. GaussPenta *gauss=NULL; GaussTria *gauss_tria=NULL; Node *node_list[numnodes]; @@ -5833,17 +5833,17 @@ /*Intermediaries */ int i,j,ig,analysis_type; - double Jdet2d,slope_magnitude,alpha2; - double xyz_list[NUMVERTICES][3]; - double xyz_list_tria[NUMVERTICES2D][3]={0.0}; - double slope[3]={0.0,0.0,0.0}; - double MAXSLOPE=.06; // 6 % - double MOUNTAINKEXPONENT=10; - double L[2][numdof]; - double DL[2][2] ={{ 0,0 },{0,0}}; //for basal drag - double DL_scalar; - double Ke_gg[numdof][numdof] ={0.0}; - double Ke_gg_gaussian[numdof][numdof]; //stiffness matrix contribution from drag + IssmDouble Jdet2d,slope_magnitude,alpha2; + IssmDouble xyz_list[NUMVERTICES][3]; + IssmDouble xyz_list_tria[NUMVERTICES2D][3]={0.0}; + IssmDouble slope[3]={0.0,0.0,0.0}; + IssmDouble MAXSLOPE=.06; // 6 % + IssmDouble MOUNTAINKEXPONENT=10; + IssmDouble L[2][numdof]; + IssmDouble DL[2][2] ={{ 0,0 },{0,0}}; //for basal drag + IssmDouble DL_scalar; + IssmDouble Ke_gg[numdof][numdof] ={0.0}; + IssmDouble Ke_gg_gaussian[numdof][numdof]; //stiffness matrix contribution from drag Friction *friction = NULL; GaussPenta *gauss=NULL; Node *node_list[numnodes]; @@ -5891,7 +5891,7 @@ slope_magnitude=sqrt(pow(slope[0],2)+pow(slope[1],2)); if (slope_magnitude>MAXSLOPE){ - alpha2=pow((double)10,MOUNTAINKEXPONENT); + alpha2=pow((IssmDouble)10,MOUNTAINKEXPONENT); } GetTriaJacobianDeterminant(&Jdet2d, &xyz_list_tria[0][0],gauss); @@ -5946,21 +5946,21 @@ /*Intermediaries */ int i,j,ig; - double Jdet; - double viscosity,stokesreconditioning; //viscosity - double epsilon[6]; /* epsilon=[exx,eyy,exy,exz,eyz];*/ - double xyz_list[NUMVERTICES][3]; - double B[4][numdofs+3]; - double Bprime[4][numdofm]; - double B2[3][numdofm]; - double Bprime2[3][numdofs+3]; - double D[4][4]={0.0}; // material matrix, simple scalar matrix. - double D2[3][3]={0.0}; // material matrix, simple scalar matrix. - double D_scalar; - double Ke_gg[numdofs][numdofm]={0.0}; //local element stiffness matrix - double Ke_gg2[numdofm][numdofs]={0.0}; //local element stiffness matrix - double Ke_gg_gaussian[numdofs+3][numdofm]; //stiffness matrix evaluated at the gaussian point. - double Ke_gg_gaussian2[numdofm][numdofs+3]; //stiffness matrix evaluated at the gaussian point. + IssmDouble Jdet; + IssmDouble viscosity,stokesreconditioning; //viscosity + IssmDouble epsilon[6]; /* epsilon=[exx,eyy,exy,exz,eyz];*/ + IssmDouble xyz_list[NUMVERTICES][3]; + IssmDouble B[4][numdofs+3]; + IssmDouble Bprime[4][numdofm]; + IssmDouble B2[3][numdofm]; + IssmDouble Bprime2[3][numdofs+3]; + IssmDouble D[4][4]={0.0}; // material matrix, simple scalar matrix. + IssmDouble D2[3][3]={0.0}; // material matrix, simple scalar matrix. + IssmDouble D_scalar; + IssmDouble Ke_gg[numdofs][numdofm]={0.0}; //local element stiffness matrix + IssmDouble Ke_gg2[numdofm][numdofs]={0.0}; //local element stiffness matrix + IssmDouble Ke_gg_gaussian[numdofs+3][numdofm]; //stiffness matrix evaluated at the gaussian point. + IssmDouble Ke_gg_gaussian2[numdofm][numdofs+3]; //stiffness matrix evaluated at the gaussian point. GaussPenta *gauss=NULL; GaussTria *gauss_tria=NULL; Node *node_list[numnodes]; @@ -6053,20 +6053,20 @@ /*Intermediaries */ int i,j,ig; int analysis_type,approximation; - double stokesreconditioning; - double viscosity,alpha2_gauss,Jdet2d; - double bed_normal[3]; - double epsilon[6]; /* epsilon=[exx,eyy,ezz,exy,exz,eyz];*/ - double xyz_list[NUMVERTICES][3]; - double xyz_list_tria[NUMVERTICES2D][3]; - double LMacAyealStokes[8][numdof2dm]; - double LprimeMacAyealStokes[8][numdof2d]; - double DLMacAyealStokes[8][8]={0.0}; - double LStokesMacAyeal[4][numdof2d]; - double LprimeStokesMacAyeal[4][numdof2dm]; - double DLStokesMacAyeal[4][4]={0.0}; - double Ke_drag_gaussian[numdof2dm][numdof2d]; - double Ke_drag_gaussian2[numdof2d][numdof2dm]; + IssmDouble stokesreconditioning; + IssmDouble viscosity,alpha2_gauss,Jdet2d; + IssmDouble bed_normal[3]; + IssmDouble epsilon[6]; /* epsilon=[exx,eyy,ezz,exy,exz,eyz];*/ + IssmDouble xyz_list[NUMVERTICES][3]; + IssmDouble xyz_list_tria[NUMVERTICES2D][3]; + IssmDouble LMacAyealStokes[8][numdof2dm]; + IssmDouble LprimeMacAyealStokes[8][numdof2d]; + IssmDouble DLMacAyealStokes[8][8]={0.0}; + IssmDouble LStokesMacAyeal[4][numdof2d]; + IssmDouble LprimeStokesMacAyeal[4][numdof2dm]; + IssmDouble DLStokesMacAyeal[4][4]={0.0}; + IssmDouble Ke_drag_gaussian[numdof2dm][numdof2d]; + IssmDouble Ke_drag_gaussian2[numdof2d][numdof2dm]; Friction* friction=NULL; GaussPenta *gauss=NULL; Node *node_list[numnodes]; @@ -6241,7 +6241,7 @@ /*Intermediaries*/ int connectivity[2]; int i,i0,i1,j0,j1; - double one0,one1; + IssmDouble one0,one1; /*Initialize Element matrix*/ ElementMatrix* Ke=new ElementMatrix(nodes,NUMVERTICES,this->parameters,NoneApproximationEnum); @@ -6258,8 +6258,8 @@ /*Find connectivity for the two nodes*/ connectivity[0]=nodes[i]->GetConnectivity(); connectivity[1]=nodes[i+3]->GetConnectivity(); - one0=1/(double)connectivity[0]; - one1=1/(double)connectivity[1]; + one0=1/(IssmDouble)connectivity[0]; + one1=1/(IssmDouble)connectivity[1]; /*Create matrix for these two nodes*/ if (IsOnBed() && IsOnSurface()){ @@ -6340,16 +6340,16 @@ /*Intermediaries */ int i,j,ig,approximation; - double Jdet; - double viscosity, oldviscosity, newviscosity, viscosity_overshoot; - double epsilon[5],oldepsilon[5]; /* epsilon=[exx,eyy,exy,exz,eyz];*/ - double epsilons[6]; //6 for stokes - double xyz_list[NUMVERTICES][3]; - double B[3][numdof2d]; - double Bprime[3][numdof2d]; - double D[3][3]={0.0}; // material matrix, simple scalar matrix. - double D_scalar; - double Ke_gg_gaussian[numdof2d][numdof2d]; //stiffness matrix evaluated at the gaussian point. + IssmDouble Jdet; + IssmDouble viscosity, oldviscosity, newviscosity, viscosity_overshoot; + IssmDouble epsilon[5],oldepsilon[5]; /* epsilon=[exx,eyy,exy,exz,eyz];*/ + IssmDouble epsilons[6]; //6 for stokes + IssmDouble xyz_list[NUMVERTICES][3]; + IssmDouble B[3][numdof2d]; + IssmDouble Bprime[3][numdof2d]; + IssmDouble D[3][3]={0.0}; // material matrix, simple scalar matrix. + IssmDouble D_scalar; + IssmDouble Ke_gg_gaussian[numdof2d][numdof2d]; //stiffness matrix evaluated at the gaussian point. Tria* tria=NULL; Penta* pentabase=NULL; GaussPenta *gauss=NULL; @@ -6493,14 +6493,14 @@ /*Intermediaries */ int i,j,ig; int approximation; - double xyz_list[NUMVERTICES][3]; - double Jdet; - double viscosity,oldviscosity,newviscosity,viscosity_overshoot; //viscosity - double epsilon[5],oldepsilon[5]; /* epsilon=[exx,eyy,exy,exz,eyz];*/ - double D_scalar; - double D[5][5]={0.0}; // material matrix, simple scalar matrix. - double B[5][numdof]; - double Bprime[5][numdof]; + IssmDouble xyz_list[NUMVERTICES][3]; + IssmDouble Jdet; + IssmDouble viscosity,oldviscosity,newviscosity,viscosity_overshoot; //viscosity + IssmDouble epsilon[5],oldepsilon[5]; /* epsilon=[exx,eyy,exy,exz,eyz];*/ + IssmDouble D_scalar; + IssmDouble D[5][5]={0.0}; // material matrix, simple scalar matrix. + IssmDouble B[5][numdof]; + IssmDouble Bprime[5][numdof]; Tria* tria=NULL; GaussPenta *gauss=NULL; @@ -6558,15 +6558,15 @@ /*Intermediaries */ int i,j,ig; int analysis_type; - double xyz_list[NUMVERTICES][3]; - double xyz_list_tria[NUMVERTICES2D][3]={0.0}; - double slope_magnitude,alpha2,Jdet; - double slope[3]={0.0,0.0,0.0}; - double MAXSLOPE=.06; // 6 % - double MOUNTAINKEXPONENT=10; - double L[2][numdof]; - double DL[2][2]={{ 0,0 },{0,0}}; //for basal drag - double DL_scalar; + IssmDouble xyz_list[NUMVERTICES][3]; + IssmDouble xyz_list_tria[NUMVERTICES2D][3]={0.0}; + IssmDouble slope_magnitude,alpha2,Jdet; + IssmDouble slope[3]={0.0,0.0,0.0}; + IssmDouble MAXSLOPE=.06; // 6 % + IssmDouble MOUNTAINKEXPONENT=10; + IssmDouble L[2][numdof]; + IssmDouble DL[2][2]={{ 0,0 },{0,0}}; //for basal drag + IssmDouble DL_scalar; Friction *friction = NULL; GaussPenta *gauss=NULL; @@ -6603,7 +6603,7 @@ // If we have a slope > 6% for this element, it means we are on a mountain. In this particular case, //velocity should be = 0. To achieve this result, we set alpha2_list to a very high value: */ if (slope_magnitude>MAXSLOPE){ - alpha2=pow((double)10,MOUNTAINKEXPONENT); + alpha2=pow((IssmDouble)10,MOUNTAINKEXPONENT); } DL_scalar=alpha2*gauss->weight*Jdet; @@ -6659,14 +6659,14 @@ /*Intermediaries */ int i,j,ig,approximation; - double Jdet,viscosity,stokesreconditioning; - double xyz_list[NUMVERTICES][3]; - double epsilon[6]; /* epsilon=[exx,eyy,ezz,exy,exz,eyz];*/ - double B[8][27]; - double B_prime[8][27]; - double D_scalar; - double D[8][8]={0.0}; - double Ke_temp[27][27]={0.0}; //for the six nodes and the bubble + IssmDouble Jdet,viscosity,stokesreconditioning; + IssmDouble xyz_list[NUMVERTICES][3]; + IssmDouble epsilon[6]; /* epsilon=[exx,eyy,ezz,exy,exz,eyz];*/ + IssmDouble B[8][27]; + IssmDouble B_prime[8][27]; + IssmDouble D_scalar; + IssmDouble D[8][8]={0.0}; + IssmDouble Ke_temp[27][27]={0.0}; //for the six nodes and the bubble GaussPenta *gauss=NULL; /*If on water or not Stokes, skip stiffness: */ @@ -6725,14 +6725,14 @@ /*Intermediaries */ int i,j,ig; int analysis_type,approximation; - double alpha2,Jdet2d; - double stokesreconditioning,viscosity; - double epsilon[6]; /* epsilon=[exx,eyy,ezz,exy,exz,eyz];*/ - double xyz_list[NUMVERTICES][3]; - double xyz_list_tria[NUMVERTICES2D][3]; - double LStokes[2][numdof2d]; - double DLStokes[2][2]={0.0}; - double Ke_drag_gaussian[numdof2d][numdof2d]; + IssmDouble alpha2,Jdet2d; + IssmDouble stokesreconditioning,viscosity; + IssmDouble epsilon[6]; /* epsilon=[exx,eyy,ezz,exy,exz,eyz];*/ + IssmDouble xyz_list[NUMVERTICES][3]; + IssmDouble xyz_list_tria[NUMVERTICES2D][3]; + IssmDouble LStokes[2][numdof2d]; + IssmDouble DLStokes[2][2]={0.0}; + IssmDouble Ke_drag_gaussian[numdof2d][numdof2d]; Friction* friction=NULL; GaussPenta *gauss=NULL; @@ -6810,11 +6810,11 @@ /*Intermediaries */ int i,j,ig; - double Jdet; - double xyz_list[NUMVERTICES][3]; - double B[NDOF1][NUMVERTICES]; - double Bprime[NDOF1][NUMVERTICES]; - double DL_scalar; + IssmDouble Jdet; + IssmDouble xyz_list[NUMVERTICES][3]; + IssmDouble B[NDOF1][NUMVERTICES]; + IssmDouble Bprime[NDOF1][NUMVERTICES]; + IssmDouble DL_scalar; GaussPenta *gauss=NULL; /*Initialize Element matrix*/ @@ -6856,11 +6856,11 @@ /*Intermediaries */ int i,j,ig; - double xyz_list[NUMVERTICES][3]; - double xyz_list_tria[NUMVERTICES2D][3]; - double surface_normal[3]; - double Jdet2d,DL_scalar; - double basis[NUMVERTICES]; + IssmDouble xyz_list[NUMVERTICES][3]; + IssmDouble xyz_list_tria[NUMVERTICES2D][3]; + IssmDouble surface_normal[3]; + IssmDouble Jdet2d,DL_scalar; + IssmDouble basis[NUMVERTICES]; GaussPenta *gauss=NULL; /*Initialize Element matrix*/ @@ -6916,13 +6916,13 @@ /*Intermediaries */ int i,j,ig; int approximation; - double viscosity,Jdet; - double stokesreconditioning; - double epsilon[6]; /* epsilon=[exx,eyy,ezz,exy,exz,eyz];*/ - double dw[3]; - double xyz_list[NUMVERTICES][3]; - double basis[6]; //for the six nodes of the penta - double dbasis[3][6]; //for the six nodes of the penta + IssmDouble viscosity,Jdet; + IssmDouble stokesreconditioning; + IssmDouble epsilon[6]; /* epsilon=[exx,eyy,ezz,exy,exz,eyz];*/ + IssmDouble dw[3]; + IssmDouble xyz_list[NUMVERTICES][3]; + IssmDouble basis[6]; //for the six nodes of the penta + IssmDouble dbasis[3][6]; //for the six nodes of the penta GaussPenta *gauss=NULL; /*Initialize Element vector and return if necessary*/ @@ -6978,15 +6978,15 @@ /*Intermediaries*/ int i,j,ig; int approximation,analysis_type; - double Jdet,Jdet2d; - double stokesreconditioning; - double bed_normal[3]; - double epsilon[6]; /* epsilon=[exx,eyy,ezz,exy,exz,eyz];*/ - double viscosity, w, alpha2_gauss; - double dw[3]; - double xyz_list_tria[NUMVERTICES2D][3]; - double xyz_list[NUMVERTICES][3]; - double basis[6]; //for the six nodes of the penta + IssmDouble Jdet,Jdet2d; + IssmDouble stokesreconditioning; + IssmDouble bed_normal[3]; + IssmDouble epsilon[6]; /* epsilon=[exx,eyy,ezz,exy,exz,eyz];*/ + IssmDouble viscosity, w, alpha2_gauss; + IssmDouble dw[3]; + IssmDouble xyz_list_tria[NUMVERTICES2D][3]; + IssmDouble xyz_list[NUMVERTICES][3]; + IssmDouble basis[6]; //for the six nodes of the penta Tria* tria=NULL; Friction* friction=NULL; GaussPenta *gauss=NULL; @@ -7067,13 +7067,13 @@ /*Intermediaries */ int i,j,ig; int approximation; - double viscosity,Jdet; - double stokesreconditioning; - double epsilon[6]; /* epsilon=[exx,eyy,ezz,exy,exz,eyz];*/ - double dw[3]; - double xyz_list[NUMVERTICES][3]; - double basis[6]; //for the six nodes of the penta - double dbasis[3][6]; //for the six nodes of the penta + IssmDouble viscosity,Jdet; + IssmDouble stokesreconditioning; + IssmDouble epsilon[6]; /* epsilon=[exx,eyy,ezz,exy,exz,eyz];*/ + IssmDouble dw[3]; + IssmDouble xyz_list[NUMVERTICES][3]; + IssmDouble basis[6]; //for the six nodes of the penta + IssmDouble dbasis[3][6]; //for the six nodes of the penta GaussPenta *gauss=NULL; /*Initialize Element vector and return if necessary*/ @@ -7129,15 +7129,15 @@ /*Intermediaries*/ int i,j,ig; int approximation,analysis_type; - double Jdet,Jdet2d; - double stokesreconditioning; - double bed_normal[3]; - double epsilon[6]; /* epsilon=[exx,eyy,ezz,exy,exz,eyz];*/ - double viscosity, w, alpha2_gauss; - double dw[3]; - double xyz_list_tria[NUMVERTICES2D][3]; - double xyz_list[NUMVERTICES][3]; - double basis[6]; //for the six nodes of the penta + IssmDouble Jdet,Jdet2d; + IssmDouble stokesreconditioning; + IssmDouble bed_normal[3]; + IssmDouble epsilon[6]; /* epsilon=[exx,eyy,ezz,exy,exz,eyz];*/ + IssmDouble viscosity, w, alpha2_gauss; + IssmDouble dw[3]; + IssmDouble xyz_list_tria[NUMVERTICES2D][3]; + IssmDouble xyz_list[NUMVERTICES][3]; + IssmDouble basis[6]; //for the six nodes of the penta Tria* tria=NULL; Friction* friction=NULL; GaussPenta *gauss=NULL; @@ -7279,14 +7279,14 @@ int i,j,k,ig; int node0,node1; int connectivity[2]; - double Jdet; - double xyz_list[NUMVERTICES][3]; - double xyz_list_segment[2][3]; - double z_list[NUMVERTICES]; - double z_segment[2],slope[2]; - double slope2,constant_part; - double rho_ice,gravity,n,B; - double ub,vb,z_g,surface,thickness; + IssmDouble Jdet; + IssmDouble xyz_list[NUMVERTICES][3]; + IssmDouble xyz_list_segment[2][3]; + IssmDouble z_list[NUMVERTICES]; + IssmDouble z_segment[2],slope[2]; + IssmDouble slope2,constant_part; + IssmDouble rho_ice,gravity,n,B; + IssmDouble ub,vb,z_g,surface,thickness; GaussPenta* gauss=NULL; /*Initialize Element vector*/ @@ -7334,22 +7334,22 @@ GetSegmentJacobianDeterminant(&Jdet,&xyz_list_segment[0][0],gauss); if (IsOnSurface()){ - for(j=0;jvalues[2*node1+j]+=constant_part*pow((surface-z_g)/B,n)*slope[j]*Jdet*gauss->weight/(double)connectivity[1]; + for(j=0;jvalues[2*node1+j]+=constant_part*pow((surface-z_g)/B,n)*slope[j]*Jdet*gauss->weight/(IssmDouble)connectivity[1]; } else{//connectivity is too large, should take only half on it - for(j=0;jvalues[2*node1+j]+=constant_part*pow((surface-z_g)/B,n)*slope[j]*Jdet*gauss->weight*2/(double)connectivity[1]; + for(j=0;jvalues[2*node1+j]+=constant_part*pow((surface-z_g)/B,n)*slope[j]*Jdet*gauss->weight*2/(IssmDouble)connectivity[1]; } } delete gauss; //Deal with lower surface if (IsOnBed()){ - constant_part=-1.58*pow((double)10.0,-(double)10.0)*rho_ice*gravity*thickness; + constant_part=-1.58*pow((IssmDouble)10.0,-(IssmDouble)10.0)*rho_ice*gravity*thickness; ub=constant_part*slope[0]; vb=constant_part*slope[1]; - pe->values[2*node0]+=ub/(double)connectivity[0]; - pe->values[2*node0+1]+=vb/(double)connectivity[0]; + pe->values[2*node0]+=ub/(IssmDouble)connectivity[0]; + pe->values[2*node0+1]+=vb/(IssmDouble)connectivity[0]; } } @@ -7379,11 +7379,11 @@ /*Intermediaries*/ int i,j,ig; - double Jdet; - double slope[3]; //do not put 2! this goes into GetInputDerivativeValue, which addresses slope[3] also! - double driving_stress_baseline,thickness; - double xyz_list[NUMVERTICES][3]; - double basis[6]; + IssmDouble Jdet; + IssmDouble slope[3]; //do not put 2! this goes into GetInputDerivativeValue, which addresses slope[3] also! + IssmDouble driving_stress_baseline,thickness; + IssmDouble xyz_list[NUMVERTICES][3]; + IssmDouble basis[6]; GaussPenta *gauss=NULL; /*Initialize Element vector*/ @@ -7442,19 +7442,19 @@ /*Intermediaries*/ int i,j,ig; int approximation; - double Jdet,viscosity; - double gravity,rho_ice,stokesreconditioning; - double xyz_list[NUMVERTICES][3]; - double epsilon[6]; /* epsilon=[exx,eyy,ezz,exy,exz,eyz];*/ - double l1l7[7]; //for the six nodes and the bubble - double B[8][numdofbubble]; - double B_prime[8][numdofbubble]; - double B_prime_bubble[8][3]; - double D[8][8]={0.0}; - double D_scalar; - double Pe_gaussian[numdofbubble]={0.0}; //for the six nodes and the bubble - double Ke_temp[numdofbubble][3]={0.0}; //for the six nodes and the bubble - double Ke_gaussian[numdofbubble][3]; + IssmDouble Jdet,viscosity; + IssmDouble gravity,rho_ice,stokesreconditioning; + IssmDouble xyz_list[NUMVERTICES][3]; + IssmDouble epsilon[6]; /* epsilon=[exx,eyy,ezz,exy,exz,eyz];*/ + IssmDouble l1l7[7]; //for the six nodes and the bubble + IssmDouble B[8][numdofbubble]; + IssmDouble B_prime[8][numdofbubble]; + IssmDouble B_prime_bubble[8][3]; + IssmDouble D[8][8]={0.0}; + IssmDouble D_scalar; + IssmDouble Pe_gaussian[numdofbubble]={0.0}; //for the six nodes and the bubble + IssmDouble Ke_temp[numdofbubble][3]={0.0}; //for the six nodes and the bubble + IssmDouble Ke_gaussian[numdofbubble][3]; GaussPenta *gauss=NULL; /*Initialize Element vector and return if necessary*/ @@ -7519,14 +7519,14 @@ /*Intermediaries*/ int i,j,ig; int approximation,shelf_dampening; - double gravity,rho_water,bed,water_pressure; - double damper,normal_vel,vx,vy,vz,dt; - double xyz_list_tria[NUMVERTICES2D][3]; - double xyz_list[NUMVERTICES][3]; - double bed_normal[3]; - double dz[3]; - double basis[6]; //for the six nodes of the penta - double Jdet2d; + IssmDouble gravity,rho_water,bed,water_pressure; + IssmDouble damper,normal_vel,vx,vy,vz,dt; + IssmDouble xyz_list_tria[NUMVERTICES2D][3]; + IssmDouble xyz_list[NUMVERTICES][3]; + IssmDouble bed_normal[3]; + IssmDouble dz[3]; + IssmDouble basis[6]; //for the six nodes of the penta + IssmDouble Jdet2d; GaussPenta *gauss=NULL; /*Initialize Element vector and return if necessary*/ @@ -7604,11 +7604,11 @@ /*Intermediaries*/ int i,ig; int approximation; - double Jdet; - double xyz_list[NUMVERTICES][3]; - double dudx,dvdy,dwdz; - double du[3],dv[3],dw[3]; - double basis[6]; + IssmDouble Jdet; + IssmDouble xyz_list[NUMVERTICES][3]; + IssmDouble dudx,dvdy,dwdz; + IssmDouble du[3],dv[3],dw[3]; + IssmDouble basis[6]; GaussPenta *gauss=NULL; /*Initialize Element vector*/ @@ -7661,12 +7661,12 @@ /*Intermediaries */ int i,j,ig; int approximation; - double xyz_list[NUMVERTICES][3]; - double xyz_list_tria[NUMVERTICES2D][3]; - double Jdet2d; - double vx,vy,vz,dbdx,dbdy,basalmeltingvalue; - double slope[3]; - double basis[NUMVERTICES]; + IssmDouble xyz_list[NUMVERTICES][3]; + IssmDouble xyz_list_tria[NUMVERTICES2D][3]; + IssmDouble Jdet2d; + IssmDouble vx,vy,vz,dbdx,dbdy,basalmeltingvalue; + IssmDouble slope[3]; + IssmDouble basis[NUMVERTICES]; GaussPenta* gauss=NULL; if (!IsOnBed()) return NULL; @@ -7767,15 +7767,15 @@ /*Intermediaries */ int i,j,ig; - double xyz_list[NUMVERTICES][3]; - double Jdet; - double eps1dotdphii,eps1dotdphij; - double eps2dotdphii,eps2dotdphij; - double mu_prime; - double epsilon[5]; /* epsilon=[exx,eyy,exy,exz,eyz];*/ - double eps1[3],eps2[3]; - double phi[NUMVERTICES]; - double dphi[3][NUMVERTICES]; + IssmDouble xyz_list[NUMVERTICES][3]; + IssmDouble Jdet; + IssmDouble eps1dotdphii,eps1dotdphij; + IssmDouble eps2dotdphii,eps2dotdphij; + IssmDouble mu_prime; + IssmDouble epsilon[5]; /* epsilon=[exx,eyy,exy,exz,eyz];*/ + IssmDouble eps1[3],eps2[3]; + IssmDouble phi[NUMVERTICES]; + IssmDouble dphi[3][NUMVERTICES]; GaussPenta *gauss=NULL; /*Initialize Jacobian with regular Pattyn (first part of the Gateau derivative)*/ @@ -7832,16 +7832,16 @@ /*Intermediaries */ int i,j,ig; - double xyz_list[NUMVERTICES][3]; - double Jdet; - double eps1dotdphii,eps1dotdphij; - double eps2dotdphii,eps2dotdphij; - double eps3dotdphii,eps3dotdphij; - double mu_prime; - double epsilon[5]; /* epsilon=[exx,eyy,exy,exz,eyz];*/ - double eps1[3],eps2[3],eps3[3]; - double phi[NUMVERTICES]; - double dphi[3][NUMVERTICES]; + IssmDouble xyz_list[NUMVERTICES][3]; + IssmDouble Jdet; + IssmDouble eps1dotdphii,eps1dotdphij; + IssmDouble eps2dotdphii,eps2dotdphij; + IssmDouble eps3dotdphii,eps3dotdphij; + IssmDouble mu_prime; + IssmDouble epsilon[5]; /* epsilon=[exx,eyy,exy,exz,eyz];*/ + IssmDouble eps1[3],eps2[3],eps3[3]; + IssmDouble phi[NUMVERTICES]; + IssmDouble dphi[3][NUMVERTICES]; GaussPenta *gauss=NULL; /*Initialize Jacobian with regular Stokes (first part of the Gateau derivative)*/ @@ -7908,8 +7908,8 @@ int i; int approximation; int* doflist=NULL; - double vx,vy; - double values[numdof]; + IssmDouble vx,vy; + IssmDouble values[numdof]; GaussPenta* gauss; /*Get approximation enum and dof list: */ @@ -7949,8 +7949,8 @@ int i; int* doflist=NULL; - double vx,vy; - double values[numdof]; + IssmDouble vx,vy; + IssmDouble values[numdof]; GaussPenta* gauss=NULL; /*Get dof list: */ @@ -7985,8 +7985,8 @@ int i; int* doflist=NULL; - double vz; - double values[numdof]; + IssmDouble vz; + IssmDouble values[numdof]; GaussPenta* gauss=NULL; /*Get dof list: */ @@ -8018,9 +8018,9 @@ int i; int* doflist=NULL; - double vx,vy,vz,p; - double stokesreconditioning; - double values[numdof]; + IssmDouble vx,vy,vz,p; + IssmDouble stokesreconditioning; + IssmDouble values[numdof]; GaussPenta *gauss; /*Get dof list: */ @@ -8057,7 +8057,7 @@ } /*}}}*/ /*FUNCTION Penta::InputUpdateFromSolutionDiagnosticHoriz {{{*/ -void Penta::InputUpdateFromSolutionDiagnosticHoriz(double* solution){ +void Penta::InputUpdateFromSolutionDiagnosticHoriz(IssmDouble* solution){ int approximation; @@ -8093,20 +8093,20 @@ } /*}}}*/ /*FUNCTION Penta::InputUpdateFromSolutionDiagnosticMacAyeal {{{*/ -void Penta::InputUpdateFromSolutionDiagnosticMacAyeal(double* solution){ +void Penta::InputUpdateFromSolutionDiagnosticMacAyeal(IssmDouble* solution){ const int numdof=NDOF2*NUMVERTICES; int i; - double rho_ice,g; - double values[numdof]; - double vx[NUMVERTICES]; - double vy[NUMVERTICES]; - double vz[NUMVERTICES]; - double vel[NUMVERTICES]; - double pressure[NUMVERTICES]; - double surface[NUMVERTICES]; - double xyz_list[NUMVERTICES][3]; + IssmDouble rho_ice,g; + IssmDouble values[numdof]; + IssmDouble vx[NUMVERTICES]; + IssmDouble vy[NUMVERTICES]; + IssmDouble vz[NUMVERTICES]; + IssmDouble vel[NUMVERTICES]; + IssmDouble pressure[NUMVERTICES]; + IssmDouble surface[NUMVERTICES]; + IssmDouble xyz_list[NUMVERTICES][3]; int *doflist = NULL; Penta *penta = NULL; @@ -8174,22 +8174,22 @@ } /*}}}*/ /*FUNCTION Penta::InputUpdateFromSolutionDiagnosticMacAyealPattyn {{{*/ -void Penta::InputUpdateFromSolutionDiagnosticMacAyealPattyn(double* solution){ +void Penta::InputUpdateFromSolutionDiagnosticMacAyealPattyn(IssmDouble* solution){ const int numdof=NDOF2*NUMVERTICES; const int numdof2d=NDOF2*NUMVERTICES2D; int i; - double rho_ice,g; - double macayeal_values[numdof]; - double pattyn_values[numdof]; - double vx[NUMVERTICES]; - double vy[NUMVERTICES]; - double vz[NUMVERTICES]; - double vel[NUMVERTICES]; - double pressure[NUMVERTICES]; - double surface[NUMVERTICES]; - double xyz_list[NUMVERTICES][3]; + IssmDouble rho_ice,g; + IssmDouble macayeal_values[numdof]; + IssmDouble pattyn_values[numdof]; + IssmDouble vx[NUMVERTICES]; + IssmDouble vy[NUMVERTICES]; + IssmDouble vz[NUMVERTICES]; + IssmDouble vel[NUMVERTICES]; + IssmDouble pressure[NUMVERTICES]; + IssmDouble surface[NUMVERTICES]; + IssmDouble xyz_list[NUMVERTICES][3]; int* doflistp = NULL; int* doflistm = NULL; Penta *penta = NULL; @@ -8258,24 +8258,24 @@ } /*}}}*/ /*FUNCTION Penta::InputUpdateFromSolutionDiagnosticMacAyealStokes {{{*/ -void Penta::InputUpdateFromSolutionDiagnosticMacAyealStokes(double* solution){ +void Penta::InputUpdateFromSolutionDiagnosticMacAyealStokes(IssmDouble* solution){ const int numdofm=NDOF2*NUMVERTICES; const int numdofs=NDOF4*NUMVERTICES; const int numdof2d=NDOF2*NUMVERTICES2D; int i; - double stokesreconditioning; - double macayeal_values[numdofm]; - double stokes_values[numdofs]; - double vx[NUMVERTICES]; - double vy[NUMVERTICES]; - double vz[NUMVERTICES]; - double vzmacayeal[NUMVERTICES]; - double vzstokes[NUMVERTICES]; - double vel[NUMVERTICES]; - double pressure[NUMVERTICES]; - double xyz_list[NUMVERTICES][3]; + IssmDouble stokesreconditioning; + IssmDouble macayeal_values[numdofm]; + IssmDouble stokes_values[numdofs]; + IssmDouble vx[NUMVERTICES]; + IssmDouble vy[NUMVERTICES]; + IssmDouble vz[NUMVERTICES]; + IssmDouble vzmacayeal[NUMVERTICES]; + IssmDouble vzstokes[NUMVERTICES]; + IssmDouble vel[NUMVERTICES]; + IssmDouble pressure[NUMVERTICES]; + IssmDouble xyz_list[NUMVERTICES][3]; int* doflistm = NULL; int* doflists = NULL; Penta *penta = NULL; @@ -8358,20 +8358,20 @@ } /*}}}*/ /*FUNCTION Penta::InputUpdateFromSolutionDiagnosticPattyn {{{*/ -void Penta::InputUpdateFromSolutionDiagnosticPattyn(double* solution){ +void Penta::InputUpdateFromSolutionDiagnosticPattyn(IssmDouble* solution){ const int numdof=NDOF2*NUMVERTICES; int i; - double rho_ice,g; - double values[numdof]; - double vx[NUMVERTICES]; - double vy[NUMVERTICES]; - double vz[NUMVERTICES]; - double vel[NUMVERTICES]; - double pressure[NUMVERTICES]; - double surface[NUMVERTICES]; - double xyz_list[NUMVERTICES][3]; + IssmDouble rho_ice,g; + IssmDouble values[numdof]; + IssmDouble vx[NUMVERTICES]; + IssmDouble vy[NUMVERTICES]; + IssmDouble vz[NUMVERTICES]; + IssmDouble vel[NUMVERTICES]; + IssmDouble pressure[NUMVERTICES]; + IssmDouble surface[NUMVERTICES]; + IssmDouble xyz_list[NUMVERTICES][3]; int* doflist = NULL; /*Get dof list: */ @@ -8432,23 +8432,23 @@ } /*}}}*/ /*FUNCTION Penta::InputUpdateFromSolutionDiagnosticPattynStokes {{{*/ -void Penta::InputUpdateFromSolutionDiagnosticPattynStokes(double* solution){ +void Penta::InputUpdateFromSolutionDiagnosticPattynStokes(IssmDouble* solution){ const int numdofp=NDOF2*NUMVERTICES; const int numdofs=NDOF4*NUMVERTICES; int i; - double pattyn_values[numdofp]; - double stokes_values[numdofs]; - double vx[NUMVERTICES]; - double vy[NUMVERTICES]; - double vz[NUMVERTICES]; - double vzpattyn[NUMVERTICES]; - double vzstokes[NUMVERTICES]; - double vel[NUMVERTICES]; - double pressure[NUMVERTICES]; - double xyz_list[NUMVERTICES][3]; - double stokesreconditioning; + IssmDouble pattyn_values[numdofp]; + IssmDouble stokes_values[numdofs]; + IssmDouble vx[NUMVERTICES]; + IssmDouble vy[NUMVERTICES]; + IssmDouble vz[NUMVERTICES]; + IssmDouble vzpattyn[NUMVERTICES]; + IssmDouble vzstokes[NUMVERTICES]; + IssmDouble vel[NUMVERTICES]; + IssmDouble pressure[NUMVERTICES]; + IssmDouble xyz_list[NUMVERTICES][3]; + IssmDouble stokesreconditioning; int* doflistp = NULL; int* doflists = NULL; Penta *penta = NULL; @@ -8526,20 +8526,20 @@ } /*}}}*/ /*FUNCTION Penta::InputUpdateFromSolutionDiagnosticHutter {{{*/ -void Penta::InputUpdateFromSolutionDiagnosticHutter(double* solution){ +void Penta::InputUpdateFromSolutionDiagnosticHutter(IssmDouble* solution){ const int numdof=NDOF2*NUMVERTICES; int i; - double rho_ice,g; - double values[numdof]; - double vx[NUMVERTICES]; - double vy[NUMVERTICES]; - double vz[NUMVERTICES]; - double vel[NUMVERTICES]; - double pressure[NUMVERTICES]; - double surface[NUMVERTICES]; - double xyz_list[NUMVERTICES][3]; + IssmDouble rho_ice,g; + IssmDouble values[numdof]; + IssmDouble vx[NUMVERTICES]; + IssmDouble vy[NUMVERTICES]; + IssmDouble vz[NUMVERTICES]; + IssmDouble vel[NUMVERTICES]; + IssmDouble pressure[NUMVERTICES]; + IssmDouble surface[NUMVERTICES]; + IssmDouble xyz_list[NUMVERTICES][3]; int* doflist = NULL; /*Get dof list: */ @@ -8589,24 +8589,24 @@ } /*}}}*/ /*FUNCTION Penta::InputUpdateFromSolutionDiagnosticVert {{{*/ -void Penta::InputUpdateFromSolutionDiagnosticVert(double* solution){ +void Penta::InputUpdateFromSolutionDiagnosticVert(IssmDouble* solution){ const int numdof=NDOF1*NUMVERTICES; int i; int approximation; - double rho_ice,g; - double values[numdof]; - double vx[NUMVERTICES]; - double vy[NUMVERTICES]; - double vz[NUMVERTICES]; - double vzmacayeal[NUMVERTICES]; - double vzpattyn[NUMVERTICES]; - double vzstokes[NUMVERTICES]; - double vel[NUMVERTICES]; - double pressure[NUMVERTICES]; - double surface[NUMVERTICES]; - double xyz_list[NUMVERTICES][3]; + IssmDouble rho_ice,g; + IssmDouble values[numdof]; + IssmDouble vx[NUMVERTICES]; + IssmDouble vy[NUMVERTICES]; + IssmDouble vz[NUMVERTICES]; + IssmDouble vzmacayeal[NUMVERTICES]; + IssmDouble vzpattyn[NUMVERTICES]; + IssmDouble vzstokes[NUMVERTICES]; + IssmDouble vel[NUMVERTICES]; + IssmDouble pressure[NUMVERTICES]; + IssmDouble surface[NUMVERTICES]; + IssmDouble xyz_list[NUMVERTICES][3]; int* doflist = NULL; @@ -8693,18 +8693,18 @@ } /*}}}*/ /*FUNCTION Penta::InputUpdateFromSolutionDiagnosticStokes {{{*/ -void Penta::InputUpdateFromSolutionDiagnosticStokes(double* solution){ +void Penta::InputUpdateFromSolutionDiagnosticStokes(IssmDouble* solution){ const int numdof=NDOF4*NUMVERTICES; int i; - double values[numdof]; - double vx[NUMVERTICES]; - double vy[NUMVERTICES]; - double vz[NUMVERTICES]; - double vel[NUMVERTICES]; - double pressure[NUMVERTICES]; - double stokesreconditioning; + IssmDouble values[numdof]; + IssmDouble vx[NUMVERTICES]; + IssmDouble vy[NUMVERTICES]; + IssmDouble vz[NUMVERTICES]; + IssmDouble vel[NUMVERTICES]; + IssmDouble pressure[NUMVERTICES]; + IssmDouble stokesreconditioning; int* doflist=NULL; /*Get dof list: */ Index: /u/astrid-r1b/morlighe/issmuci/trunk-jpl/../trunk-jpl/src/c/objects/Elements/PentaRef.h =================================================================== --- /u/astrid-r1b/morlighe/issmuci/trunk-jpl/../trunk-jpl/src/c/objects/Elements/PentaRef.h (revision 12470) +++ /u/astrid-r1b/morlighe/issmuci/trunk-jpl/../trunk-jpl/src/c/objects/Elements/PentaRef.h (revision 12471) @@ -22,42 +22,42 @@ void SetElementType(int type,int type_counter); /*Numerics*/ - void GetNodalFunctionsP1(double* l1l6, GaussPenta* gauss); - void GetNodalFunctionsMINI(double* l1l7, GaussPenta* gauss); - void GetNodalFunctionsP1Derivatives(double* dh1dh6,double* xyz_list, GaussPenta* gauss); - void GetNodalFunctionsMINIDerivatives(double* dh1dh7,double* xyz_list, GaussPenta* gauss); - void GetNodalFunctionsP1DerivativesReference(double* dl1dl6,GaussPenta* gauss); - void GetNodalFunctionsMINIDerivativesReference(double* dl1dl7,GaussPenta* gauss); - void GetQuadNodalFunctions(double* l1l4,GaussPenta* gauss,int index1,int index2,int index3,int index4); - void GetQuadJacobianDeterminant(double* Jdet, double xyz_list[4][3],GaussPenta* gauss); - void GetJacobian(double* J, double* xyz_list,GaussPenta* gauss); - void GetJacobianDeterminant(double* Jdet, double* xyz_list,GaussPenta* gauss); - void GetTriaJacobianDeterminant(double* Jdet, double* xyz_list,GaussPenta* gauss); - void GetSegmentJacobianDeterminant(double* Jdet, double* xyz_list,GaussPenta* gauss); - void GetJacobianInvert(double* Jinv, double* xyz_list,GaussPenta* gauss); - void GetBMacAyealPattyn(double* B, double* xyz_list, GaussPenta* gauss); - void GetBMacAyealStokes(double* B, double* xyz_list, GaussPenta* gauss); - void GetBPattyn(double* B, double* xyz_list, GaussPenta* gauss); - void GetBStokes(double* B, double* xyz_list, GaussPenta* gauss); - void GetBprimeMacAyealStokes(double* Bprime, double* xyz_list, GaussPenta* gauss); - void GetBprimePattyn(double* B, double* xyz_list, GaussPenta* gauss); - void GetBprimeStokes(double* B_prime, double* xyz_list, GaussPenta* gauss); - void GetBprimeVert(double* B, double* xyz_list, GaussPenta* gauss); - void GetBAdvec(double* B_advec, double* xyz_list, GaussPenta* gauss); - void GetBConduct(double* B_conduct, double* xyz_list, GaussPenta* gauss); - void GetBVert(double* B, double* xyz_list, GaussPenta* gauss); - void GetBprimeAdvec(double* Bprime_advec, double* xyz_list, GaussPenta* gauss); - void GetL(double* L, GaussPenta* gauss,int numdof); - void GetLStokes(double* LStokes, GaussPenta* gauss); - void GetLprimeStokes(double* LprimeStokes, double* xyz_list, GaussPenta* gauss); - void GetLMacAyealStokes(double* LMacAyealStokes, GaussPenta* gauss); - void GetLprimeMacAyealStokes(double* LprimeMacAyealStokes, double* xyz_list, GaussPenta* gauss); - void GetLStokesMacAyeal(double* LStokesMacAyeal, GaussPenta* gauss); - void GetLprimeStokesMacAyeal(double* LprimeStokesMacAyeal, double* xyz_list, GaussPenta* gauss); - void GetInputValue(double* pvalue,double* plist, GaussPenta* gauss); - void GetInputValue(double* pvalue,double* plist,GaussTria* gauss){_error_("only PentaGauss are supported");}; - void GetInputDerivativeValue(double* pvalues, double* plist,double* xyz_list, GaussPenta* gauss); - void GetInputDerivativeValue(double* pvalues, double* plist,double* xyz_list, GaussTria* gauss){_error_("only PentaGauss are supported");}; + void GetNodalFunctionsP1(IssmDouble* l1l6, GaussPenta* gauss); + void GetNodalFunctionsMINI(IssmDouble* l1l7, GaussPenta* gauss); + void GetNodalFunctionsP1Derivatives(IssmDouble* dh1dh6,IssmDouble* xyz_list, GaussPenta* gauss); + void GetNodalFunctionsMINIDerivatives(IssmDouble* dh1dh7,IssmDouble* xyz_list, GaussPenta* gauss); + void GetNodalFunctionsP1DerivativesReference(IssmDouble* dl1dl6,GaussPenta* gauss); + void GetNodalFunctionsMINIDerivativesReference(IssmDouble* dl1dl7,GaussPenta* gauss); + void GetQuadNodalFunctions(IssmDouble* l1l4,GaussPenta* gauss,int index1,int index2,int index3,int index4); + void GetQuadJacobianDeterminant(IssmDouble* Jdet, IssmDouble xyz_list[4][3],GaussPenta* gauss); + void GetJacobian(IssmDouble* J, IssmDouble* xyz_list,GaussPenta* gauss); + void GetJacobianDeterminant(IssmDouble* Jdet, IssmDouble* xyz_list,GaussPenta* gauss); + void GetTriaJacobianDeterminant(IssmDouble* Jdet, IssmDouble* xyz_list,GaussPenta* gauss); + void GetSegmentJacobianDeterminant(IssmDouble* Jdet, IssmDouble* xyz_list,GaussPenta* gauss); + void GetJacobianInvert(IssmDouble* Jinv, IssmDouble* xyz_list,GaussPenta* gauss); + void GetBMacAyealPattyn(IssmDouble* B, IssmDouble* xyz_list, GaussPenta* gauss); + void GetBMacAyealStokes(IssmDouble* B, IssmDouble* xyz_list, GaussPenta* gauss); + void GetBPattyn(IssmDouble* B, IssmDouble* xyz_list, GaussPenta* gauss); + void GetBStokes(IssmDouble* B, IssmDouble* xyz_list, GaussPenta* gauss); + void GetBprimeMacAyealStokes(IssmDouble* Bprime, IssmDouble* xyz_list, GaussPenta* gauss); + void GetBprimePattyn(IssmDouble* B, IssmDouble* xyz_list, GaussPenta* gauss); + void GetBprimeStokes(IssmDouble* B_prime, IssmDouble* xyz_list, GaussPenta* gauss); + void GetBprimeVert(IssmDouble* B, IssmDouble* xyz_list, GaussPenta* gauss); + void GetBAdvec(IssmDouble* B_advec, IssmDouble* xyz_list, GaussPenta* gauss); + void GetBConduct(IssmDouble* B_conduct, IssmDouble* xyz_list, GaussPenta* gauss); + void GetBVert(IssmDouble* B, IssmDouble* xyz_list, GaussPenta* gauss); + void GetBprimeAdvec(IssmDouble* Bprime_advec, IssmDouble* xyz_list, GaussPenta* gauss); + void GetL(IssmDouble* L, GaussPenta* gauss,int numdof); + void GetLStokes(IssmDouble* LStokes, GaussPenta* gauss); + void GetLprimeStokes(IssmDouble* LprimeStokes, IssmDouble* xyz_list, GaussPenta* gauss); + void GetLMacAyealStokes(IssmDouble* LMacAyealStokes, GaussPenta* gauss); + void GetLprimeMacAyealStokes(IssmDouble* LprimeMacAyealStokes, IssmDouble* xyz_list, GaussPenta* gauss); + void GetLStokesMacAyeal(IssmDouble* LStokesMacAyeal, GaussPenta* gauss); + void GetLprimeStokesMacAyeal(IssmDouble* LprimeStokesMacAyeal, IssmDouble* xyz_list, GaussPenta* gauss); + void GetInputValue(IssmDouble* pvalue,IssmDouble* plist, GaussPenta* gauss); + void GetInputValue(IssmDouble* pvalue,IssmDouble* plist,GaussTria* gauss){_error_("only PentaGauss are supported");}; + void GetInputDerivativeValue(IssmDouble* pvalues, IssmDouble* plist,IssmDouble* xyz_list, GaussPenta* gauss); + void GetInputDerivativeValue(IssmDouble* pvalues, IssmDouble* plist,IssmDouble* xyz_list, GaussTria* gauss){_error_("only PentaGauss are supported");}; }; #endif Index: /u/astrid-r1b/morlighe/issmuci/trunk-jpl/../trunk-jpl/src/c/objects/Elements/Penta.h =================================================================== --- /u/astrid-r1b/morlighe/issmuci/trunk-jpl/../trunk-jpl/src/c/objects/Elements/Penta.h (revision 12470) +++ /u/astrid-r1b/morlighe/issmuci/trunk-jpl/../trunk-jpl/src/c/objects/Elements/Penta.h (revision 12471) @@ -58,22 +58,22 @@ /*}}}*/ /*Update virtual functions definitions: {{{*/ void InputUpdateFromConstant(bool constant, int name); - void InputUpdateFromConstant(double constant, int name); + void InputUpdateFromConstant(IssmDouble constant, int name); void InputUpdateFromConstant(int constant, int name); - void InputUpdateFromSolution(double* solutiong); + void InputUpdateFromSolution(IssmDouble* solutiong); void InputUpdateFromVector(bool* vector, int name, int type); - void InputUpdateFromVector(double* vector, int name, int type); + void InputUpdateFromVector(IssmDouble* vector, int name, int type); void InputUpdateFromVector(int* vector, int name, int type); #ifdef _HAVE_DAKOTA_ void InputUpdateFromVectorDakota(bool* vector, int name, int type); - void InputUpdateFromVectorDakota(double* vector, int name, int type); + void InputUpdateFromVectorDakota(IssmDouble* vector, int name, int type); void InputUpdateFromVectorDakota(int* vector, int name, int type); - void InputUpdateFromMatrixDakota(double* matrix, int nows, int ncols, int name, int type); + void InputUpdateFromMatrixDakota(IssmDouble* matrix, int nows, int ncols, int name, int type); #endif void InputUpdateFromIoModel(int index, IoModel* iomodel); /*}}}*/ /*Element virtual functions definitions: {{{*/ - void AverageOntoPartition(Vector* partition_contributions,Vector* partition_areas,double* vertex_response,double* qmu_part); + void AverageOntoPartition(Vector* partition_contributions,Vector* partition_areas,IssmDouble* vertex_response,IssmDouble* qmu_part); void BasalFrictionCreateInput(void); void ComputeBasalStress(Vector* sigma_b); void ComputeStrainRate(Vector* eps); @@ -86,57 +86,57 @@ void DeleteResults(void); int GetNodeIndex(Node* node); void GetSolutionFromInputs(Vector* solution); - double GetZcoord(GaussPenta* gauss); + IssmDouble GetZcoord(GaussPenta* gauss); void GetVectorFromInputs(Vector* vector,int name_enum); void GetVectorFromResults(Vector* vector,int offset,int interp); int Sid(); - void InputArtificialNoise(int enum_type,double min, double max); - bool InputConvergence(double* eps, int* enums,int num_enums,int* criterionenums,double* criterionvalues,int num_criterionenums); - void InputCreate(double scalar,int name,int code); - void InputCreate(double* vector, int index,IoModel* iomodel,int M,int N,int vector_type,int vector_enum,int code); + void InputArtificialNoise(int enum_type,IssmDouble min, IssmDouble max); + bool InputConvergence(IssmDouble* eps, int* enums,int num_enums,int* criterionenums,IssmDouble* criterionvalues,int num_criterionenums); + void InputCreate(IssmDouble scalar,int name,int code); + void InputCreate(IssmDouble* vector, int index,IoModel* iomodel,int M,int N,int vector_type,int vector_enum,int code); void InputDepthAverageAtBase(int enum_type,int average_enum_type,int object_enum=MeshElementsEnum); void InputDuplicate(int original_enum,int new_enum); - void InputScale(int enum_type,double scale_factor); + void InputScale(int enum_type,IssmDouble scale_factor); - void InputToResult(int enum_type,int step,double time); - void MigrateGroundingLine(double* old_floating_ice,double* sheet_ungrounding); + void InputToResult(int enum_type,int step,IssmDouble time); + void MigrateGroundingLine(IssmDouble* old_floating_ice,IssmDouble* sheet_ungrounding); void PotentialSheetUngrounding(Vector* potential_sheet_ungrounding); - void RequestedOutput(int output_enum,int step,double time); - void ListResultsInfo(int** results_enums,int** results_size,double** results_times,int** results_steps,int* num_results); + void RequestedOutput(int output_enum,int step,IssmDouble time); + void ListResultsInfo(int** results_enums,int** results_size,IssmDouble** results_times,int** results_steps,int* num_results); void PatchFill(int* pcount, Patch* patch); void PatchSize(int* pnumrows, int* pnumvertices,int* pnumnodes); - void PositiveDegreeDay(double* pdds,double* pds,double signorm); + void PositiveDegreeDay(IssmDouble* pdds,IssmDouble* pds,IssmDouble signorm); void ProcessResultsUnits(void); void ResetCoordinateSystem(void); - double SurfaceArea(void); + IssmDouble SurfaceArea(void); void Update(int index, IoModel* iomodel,int analysis_counter,int analysis_type); - int UpdatePotentialSheetUngrounding(double* potential_sheet_ungrounding,Vector* vec_nodes_on_iceshelf,double* nodes_on_iceshelf); - int NodalValue(double* pvalue, int index, int natureofdataenum,bool process_units); - double TimeAdapt(); + int UpdatePotentialSheetUngrounding(IssmDouble* potential_sheet_ungrounding,Vector* vec_nodes_on_iceshelf,IssmDouble* nodes_on_iceshelf); + int NodalValue(IssmDouble* pvalue, int index, int natureofdataenum,bool process_units); + IssmDouble TimeAdapt(); int* GetHorizontalNeighboorSids(void); void ViscousHeatingCreateInput(void); - void SmearFunction(Vector* smearedvector,double (*WeightFunction)(double distance,double radius),double radius); + void SmearFunction(Vector* smearedvector,IssmDouble (*WeightFunction)(IssmDouble distance,IssmDouble radius),IssmDouble radius); #ifdef _HAVE_RESPONSES_ - double IceVolume(void); - 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 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 ElementResponse(double* presponse,int response_enum,bool process_units); - void MaxVx(double* pmaxvx, bool process_units); - void MaxVy(double* pmaxvy, bool process_units); - void MaxVz(double* pmaxvz, bool process_units); + IssmDouble IceVolume(void); + void MinVel(IssmDouble* pminvel, bool process_units); + void MinVx(IssmDouble* pminvx, bool process_units); + void MinVy(IssmDouble* pminvy, bool process_units); + void MinVz(IssmDouble* pminvz, bool process_units); + IssmDouble MassFlux(IssmDouble* segment,bool process_units); + void MaxAbsVx(IssmDouble* pmaxabsvx, bool process_units); + void MaxAbsVy(IssmDouble* pmaxabsvy, bool process_units); + void MaxAbsVz(IssmDouble* pmaxabsvz, bool process_units); + void MaxVel(IssmDouble* pmaxvel, bool process_units); + void ElementResponse(IssmDouble* presponse,int response_enum,bool process_units); + void MaxVx(IssmDouble* pmaxvx, bool process_units); + void MaxVy(IssmDouble* pmaxvy, bool process_units); + void MaxVz(IssmDouble* pmaxvz, bool process_units); #endif #ifdef _HAVE_CONTROL_ - double DragCoefficientAbsGradient(bool process_units,int weight_index); + IssmDouble DragCoefficientAbsGradient(bool process_units,int weight_index); void GradientIndexing(int* indexing,int control_index); void Gradj(Vector* gradient,int control_type,int control_index); void GradjDragMacAyeal(Vector* gradient,int control_index); @@ -146,23 +146,23 @@ void GradjBbarPattyn(Vector* gradient,int control_index); void GradjBbarStokes(Vector* gradient,int control_index); void GetVectorFromControlInputs(Vector* gradient,int control_enum,int control_index,const char* data); - void SetControlInputsFromVector(double* vector,int control_enum,int control_index); + void SetControlInputsFromVector(IssmDouble* vector,int control_enum,int control_index); void ControlInputGetGradient(Vector* gradient,int enum_type,int control_index); - void ControlInputScaleGradient(int enum_type,double scale); - void ControlInputSetGradient(double* gradient,int enum_type,int control_index); - double RheologyBbarAbsGradient(bool process_units,int weight_index); - double ThicknessAbsMisfit( bool process_units,int weight_index); - double SurfaceAbsVelMisfit( bool process_units,int weight_index); - double SurfaceRelVelMisfit( bool process_units,int weight_index); - double SurfaceLogVelMisfit( bool process_units,int weight_index); - double SurfaceLogVxVyMisfit( bool process_units,int weight_index); - double SurfaceAverageVelMisfit(bool process_units,int weight_index); - double ThicknessAbsGradient(bool process_units,int weight_index); - void InputControlUpdate(double scalar,bool save_parameter); + void ControlInputScaleGradient(int enum_type,IssmDouble scale); + void ControlInputSetGradient(IssmDouble* gradient,int enum_type,int control_index); + IssmDouble RheologyBbarAbsGradient(bool process_units,int weight_index); + IssmDouble ThicknessAbsMisfit( bool process_units,int weight_index); + IssmDouble SurfaceAbsVelMisfit( bool process_units,int weight_index); + IssmDouble SurfaceRelVelMisfit( bool process_units,int weight_index); + IssmDouble SurfaceLogVelMisfit( bool process_units,int weight_index); + IssmDouble SurfaceLogVxVyMisfit( bool process_units,int weight_index); + IssmDouble SurfaceAverageVelMisfit(bool process_units,int weight_index); + IssmDouble ThicknessAbsGradient(bool process_units,int weight_index); + void InputControlUpdate(IssmDouble scalar,bool save_parameter); #endif /*}}}*/ /*Penta specific routines:{{{*/ - void BedNormal(double* bed_normal, double xyz_list[3][3]); + void BedNormal(IssmDouble* bed_normal, IssmDouble xyz_list[3][3]); ElementMatrix* CreateKMatrixPrognostic(void); ElementMatrix* CreateKMatrixSlope(void); ElementVector* CreatePVectorPrognostic(void); @@ -172,35 +172,35 @@ void GetSidList(int* sidlist); void GetConnectivityList(int* connectivity); int GetElementType(void); - void GetElementSizes(double* hx,double* hy,double* hz); - void GetInputListOnVertices(double* pvalue,int enumtype); - void GetInputListOnVertices(double* pvalue,int enumtype,double defaultvalue); - void GetInputValue(double* pvalue,Node* node,int enumtype); - void GetPhi(double* phi, double* epsilon, double viscosity); + void GetElementSizes(IssmDouble* hx,IssmDouble* hy,IssmDouble* hz); + void GetInputListOnVertices(IssmDouble* pvalue,int enumtype); + void GetInputListOnVertices(IssmDouble* pvalue,int enumtype,IssmDouble defaultvalue); + void GetInputValue(IssmDouble* pvalue,Node* node,int enumtype); + void GetPhi(IssmDouble* phi, IssmDouble* epsilon, IssmDouble viscosity); void GetSolutionFromInputsEnthalpy(Vector* solutiong); - double GetStabilizationParameter(double u, double v, double w, double diameter, double kappa); - void GetStrainRate3dPattyn(double* epsilon,double* xyz_list, GaussPenta* gauss, Input* vx_input, Input* vy_input); - void GetStrainRate3d(double* epsilon,double* xyz_list, GaussPenta* gauss, Input* vx_input, Input* vy_input, Input* vz_input); + IssmDouble GetStabilizationParameter(IssmDouble u, IssmDouble v, IssmDouble w, IssmDouble diameter, IssmDouble kappa); + void GetStrainRate3dPattyn(IssmDouble* epsilon,IssmDouble* xyz_list, GaussPenta* gauss, Input* vx_input, Input* vy_input); + void GetStrainRate3d(IssmDouble* epsilon,IssmDouble* xyz_list, GaussPenta* gauss, Input* vx_input, Input* vy_input, Input* vz_input); Penta* GetUpperElement(void); Penta* GetLowerElement(void); Penta* GetBasalElement(void); void InputExtrude(int enum_type,int object_type); - void InputUpdateFromSolutionPrognostic(double* solutiong); - void InputUpdateFromSolutionOneDof(double* solutiong,int enum_type); - void InputUpdateFromSolutionOneDofCollapsed(double* solutiong,int enum_type); + void InputUpdateFromSolutionPrognostic(IssmDouble* solutiong); + void InputUpdateFromSolutionOneDof(IssmDouble* solutiong,int enum_type); + void InputUpdateFromSolutionOneDofCollapsed(IssmDouble* solutiong,int enum_type); bool IsInput(int name); bool IsOnSurface(void); bool IsOnBed(void); bool IsFloating(void); bool IsNodeOnShelf(); - bool IsNodeOnShelfFromFlags(double* flags); + bool IsNodeOnShelfFromFlags(IssmDouble* flags); bool IsOnWater(void); - double MinEdgeLength(double xyz_list[6][3]); - void ReduceMatrixStokes(double* Ke_reduced, double* Ke_temp); - void ReduceVectorStokes(double* Pe_reduced, double* Ke_temp, double* Pe_temp); + IssmDouble MinEdgeLength(IssmDouble xyz_list[6][3]); + void ReduceMatrixStokes(IssmDouble* Ke_reduced, IssmDouble* Ke_temp); + void ReduceVectorStokes(IssmDouble* Pe_reduced, IssmDouble* Ke_temp, IssmDouble* Pe_temp); void SetClone(int* minranks); Tria* SpawnTria(int g0, int g1, int g2); - void SurfaceNormal(double* surface_normal, double xyz_list[3][3]); + void SurfaceNormal(IssmDouble* surface_normal, IssmDouble xyz_list[3][3]); #ifdef _HAVE_DIAGNOSTIC_ ElementMatrix* CreateKMatrixCouplingMacAyealPattyn(void); @@ -235,15 +235,15 @@ ElementMatrix* CreateJacobianDiagnosticMacayeal2d(void); ElementMatrix* CreateJacobianDiagnosticPattyn(void); ElementMatrix* CreateJacobianDiagnosticStokes(void); - void InputUpdateFromSolutionDiagnosticHoriz( double* solutiong); - void InputUpdateFromSolutionDiagnosticMacAyeal( double* solutiong); - void InputUpdateFromSolutionDiagnosticMacAyealPattyn( double* solutiong); - void InputUpdateFromSolutionDiagnosticMacAyealStokes( double* solutiong); - void InputUpdateFromSolutionDiagnosticPattyn( double* solutiong); - void InputUpdateFromSolutionDiagnosticPattynStokes( double* solutiong); - void InputUpdateFromSolutionDiagnosticHutter( double* solutiong); - void InputUpdateFromSolutionDiagnosticVert( double* solutiong); - void InputUpdateFromSolutionDiagnosticStokes( double* solutiong); + void InputUpdateFromSolutionDiagnosticHoriz( IssmDouble* solutiong); + void InputUpdateFromSolutionDiagnosticMacAyeal( IssmDouble* solutiong); + void InputUpdateFromSolutionDiagnosticMacAyealPattyn( IssmDouble* solutiong); + void InputUpdateFromSolutionDiagnosticMacAyealStokes( IssmDouble* solutiong); + void InputUpdateFromSolutionDiagnosticPattyn( IssmDouble* solutiong); + void InputUpdateFromSolutionDiagnosticPattynStokes( IssmDouble* solutiong); + void InputUpdateFromSolutionDiagnosticHutter( IssmDouble* solutiong); + void InputUpdateFromSolutionDiagnosticVert( IssmDouble* solutiong); + void InputUpdateFromSolutionDiagnosticStokes( IssmDouble* solutiong); void GetSolutionFromInputsDiagnosticHoriz(Vector* solutiong); void GetSolutionFromInputsDiagnosticHutter(Vector* solutiong); void GetSolutionFromInputsDiagnosticStokes(Vector* solutiong); @@ -277,8 +277,8 @@ ElementVector* CreatePVectorAdjointMacAyeal(void); ElementVector* CreatePVectorAdjointPattyn(void); ElementVector* CreatePVectorAdjointStokes(void); - void InputUpdateFromSolutionAdjointHoriz( double* solutiong); - void InputUpdateFromSolutionAdjointStokes( double* solutiong); + void InputUpdateFromSolutionAdjointHoriz( IssmDouble* solutiong); + void InputUpdateFromSolutionAdjointStokes( IssmDouble* solutiong); #endif #ifdef _HAVE_HYDROLOGY_ @@ -302,8 +302,8 @@ ElementVector* CreatePVectorThermalShelf(void); ElementVector* CreatePVectorThermalSheet(void); void GetSolutionFromInputsThermal(Vector* solutiong); - void InputUpdateFromSolutionThermal( double* solutiong); - void InputUpdateFromSolutionEnthalpy( double* solutiong); + void InputUpdateFromSolutionThermal( IssmDouble* solutiong); + void InputUpdateFromSolutionEnthalpy( IssmDouble* solutiong); #endif #ifdef _HAVE_BALANCED_ ElementMatrix* CreateKMatrixBalancethickness(void); Index: /u/astrid-r1b/morlighe/issmuci/trunk-jpl/../trunk-jpl/src/c/objects/Elements/TriaRef.cpp =================================================================== --- /u/astrid-r1b/morlighe/issmuci/trunk-jpl/../trunk-jpl/src/c/objects/Elements/TriaRef.cpp (revision 12470) +++ /u/astrid-r1b/morlighe/issmuci/trunk-jpl/../trunk-jpl/src/c/objects/Elements/TriaRef.cpp (revision 12471) @@ -56,7 +56,7 @@ /*Reference Element numerics*/ /*FUNCTION TriaRef::GetBMacAyeal {{{*/ -void TriaRef::GetBMacAyeal(double* B, double* xyz_list, GaussTria* gauss){ +void TriaRef::GetBMacAyeal(IssmDouble* B, IssmDouble* xyz_list, GaussTria* gauss){ /*Compute B matrix. B=[B1 B2 B3] where Bi is of size 3*NDOF2. * For node i, Bi can be expressed in the actual coordinate system * by: @@ -69,7 +69,7 @@ */ int i; - double dbasis[NDOF2][NUMNODES]; + IssmDouble dbasis[NDOF2][NUMNODES]; /*Get dh1dh2dh3 in actual coordinate system: */ GetNodalFunctionsDerivatives(&dbasis[0][0],xyz_list,gauss); @@ -86,7 +86,7 @@ } /*}}}*/ /*FUNCTION TriaRef::GetBMacAyealStokes {{{*/ -void TriaRef::GetBMacAyealStokes(double* B, double* xyz_list, GaussTria* gauss){ +void TriaRef::GetBMacAyealStokes(IssmDouble* B, IssmDouble* xyz_list, GaussTria* gauss){ /*Compute B matrix. B=[B1 B2 B3] where Bi is of size 3*NDOF2. * For node i, Bi can be expressed in the actual coordinate system @@ -100,7 +100,7 @@ */ /*Same thing in the actual coordinate system: */ - double dbasis[NDOF2][NUMNODES]; + IssmDouble dbasis[NDOF2][NUMNODES]; /*Get dh1dh2dh3 in actual coordinates system : */ GetNodalFunctionsDerivatives(&dbasis[0][0],xyz_list,gauss); @@ -117,7 +117,7 @@ } /*}}}*/ /*FUNCTION TriaRef::GetSegmentBFlux{{{*/ -void TriaRef::GetSegmentBFlux(double* B,GaussTria* gauss, int index1,int index2){ +void TriaRef::GetSegmentBFlux(IssmDouble* B,GaussTria* gauss, int index1,int index2){ /*Compute B matrix. B=[phi1 phi2 -phi3 -phi4] * * and phi1=phi3 phi2=phi4 @@ -125,7 +125,7 @@ * We assume B has been allocated already, of size: 1x4 */ - double l1l3[NUMNODES]; + IssmDouble l1l3[NUMNODES]; GetNodalFunctions(&l1l3[0],gauss); @@ -136,7 +136,7 @@ } /*}}}*/ /*FUNCTION TriaRef::GetSegmentBprimeFlux{{{*/ -void TriaRef::GetSegmentBprimeFlux(double* Bprime,GaussTria* gauss, int index1,int index2){ +void TriaRef::GetSegmentBprimeFlux(IssmDouble* Bprime,GaussTria* gauss, int index1,int index2){ /*Compute Bprime matrix. Bprime=[phi1 phi2 phi3 phi4] * * and phi1=phi3 phi2=phi4 @@ -144,7 +144,7 @@ * We assume Bprime has been allocated already, of size: 1x4 */ - double l1l3[NUMNODES]; + IssmDouble l1l3[NUMNODES]; GetNodalFunctions(&l1l3[0],gauss); @@ -155,7 +155,7 @@ } /*}}}*/ /*FUNCTION TriaRef::GetBPrognostic{{{*/ -void TriaRef::GetBPrognostic(double* B_prog, double* xyz_list, GaussTria* gauss){ +void TriaRef::GetBPrognostic(IssmDouble* B_prog, IssmDouble* xyz_list, GaussTria* gauss){ /*Compute B matrix. B=[B1 B2 B3] where Bi is of size 3*NDOF2. * For node i, Bi can be expressed in the actual coordinate system * by: @@ -166,7 +166,7 @@ * We assume B_prog has been allocated already, of size: 2x(NDOF1*NUMNODES) */ - double basis[NUMNODES]; + IssmDouble basis[NUMNODES]; /*Get dh1dh2dh3 in actual coordinate system: */ GetNodalFunctions(&basis[0],gauss); @@ -179,7 +179,7 @@ } /*}}}*/ /*FUNCTION TriaRef::GetBprimeMacAyeal {{{*/ -void TriaRef::GetBprimeMacAyeal(double* Bprime, double* xyz_list, GaussTria* gauss){ +void TriaRef::GetBprimeMacAyeal(IssmDouble* Bprime, IssmDouble* xyz_list, GaussTria* gauss){ /*Compute B' matrix. B'=[B1' B2' B3'] where Bi' is of size 3*NDOF2. * For node i, Bi' can be expressed in the actual coordinate system @@ -193,7 +193,7 @@ */ /*Same thing in the actual coordinate system: */ - double dbasis[NDOF2][NUMNODES]; + IssmDouble dbasis[NDOF2][NUMNODES]; /*Get dh1dh2dh3 in actual coordinates system : */ GetNodalFunctionsDerivatives(&dbasis[0][0],xyz_list,gauss); @@ -210,7 +210,7 @@ } /*}}}*/ /*FUNCTION TriaRef::GetBprimeMacAyealStokes {{{*/ -void TriaRef::GetBprimeMacAyealStokes(double* Bprime, double* xyz_list, GaussTria* gauss){ +void TriaRef::GetBprimeMacAyealStokes(IssmDouble* Bprime, IssmDouble* xyz_list, GaussTria* gauss){ /*Compute Bprime matrix. Bprime=[Bprime1 Bprime2 Bprime3] where Bprimei is of size 3*NDOF2. * For node i, Bprimei can be expressed in the actual coordinate system @@ -225,7 +225,7 @@ */ /*Same thing in the actual coordinate system: */ - double dbasis[NDOF2][NUMNODES]; + IssmDouble dbasis[NDOF2][NUMNODES]; /*Get dh1dh2dh3 in actual coordinates system : */ GetNodalFunctionsDerivatives(&dbasis[0][0],xyz_list,gauss); @@ -244,7 +244,7 @@ } /*}}}*/ /*FUNCTION TriaRef::GetBprimePrognostic{{{*/ -void TriaRef::GetBprimePrognostic(double* Bprime_prog, double* xyz_list, GaussTria* gauss){ +void TriaRef::GetBprimePrognostic(IssmDouble* Bprime_prog, IssmDouble* xyz_list, GaussTria* gauss){ /*Compute B' matrix. B'=[B1' B2' B3'] where Bi' is of size 3*NDOF2. * For node i, Bi' can be expressed in the actual coordinate system * by: @@ -256,7 +256,7 @@ */ /*Same thing in the actual coordinate system: */ - double dbasis[NDOF2][NUMNODES]; + IssmDouble dbasis[NDOF2][NUMNODES]; /*Get dh1dh2dh3 in actual coordinates system : */ GetNodalFunctionsDerivatives(&dbasis[0][0],xyz_list,gauss); @@ -269,7 +269,7 @@ } /*}}}*/ /*FUNCTION TriaRef::GetL{{{*/ -void TriaRef::GetL(double* L, double* xyz_list,GaussTria* gauss,int numdof){ +void TriaRef::GetL(IssmDouble* L, IssmDouble* xyz_list,GaussTria* gauss,int numdof){ /*Compute L matrix. L=[L1 L2 L3] where Li is square and of size numdof. * For node i, Li can be expressed in the actual coordinate system * by: @@ -284,7 +284,7 @@ */ int i; - double basis[3]; + IssmDouble basis[3]; /*Get basis in actual coordinate system: */ GetNodalFunctions(basis,gauss); @@ -306,10 +306,10 @@ } /*}}}*/ /*FUNCTION TriaRef::GetJacobian{{{*/ -void TriaRef::GetJacobian(double* J, double* xyz_list,GaussTria* gauss){ +void TriaRef::GetJacobian(IssmDouble* J, IssmDouble* xyz_list,GaussTria* gauss){ /*The Jacobian is constant over the element, discard the gaussian points. * J is assumed to have been allocated of size NDOF2xNDOF2.*/ - double x1,y1,x2,y2,x3,y3; + IssmDouble x1,y1,x2,y2,x3,y3; x1=*(xyz_list+NUMNODES*0+0); y1=*(xyz_list+NUMNODES*0+1); @@ -326,10 +326,10 @@ } /*}}}*/ /*FUNCTION TriaRef::GetSegmentJacobianDeterminant{{{*/ -void TriaRef::GetSegmentJacobianDeterminant(double* Jdet, double* xyz_list,GaussTria* gauss){ +void TriaRef::GetSegmentJacobianDeterminant(IssmDouble* Jdet, IssmDouble* xyz_list,GaussTria* gauss){ /*The Jacobian determinant is constant over the element, discard the gaussian points. * J is assumed to have been allocated*/ - double x1,y1,x2,y2; + IssmDouble x1,y1,x2,y2; x1=*(xyz_list+3*0+0); y1=*(xyz_list+3*0+1); @@ -342,10 +342,10 @@ } /*}}}*/ /*FUNCTION TriaRef::GetJacobianDeterminant2d{{{*/ -void TriaRef::GetJacobianDeterminant2d(double* Jdet, double* xyz_list,GaussTria* gauss){ +void TriaRef::GetJacobianDeterminant2d(IssmDouble* Jdet, IssmDouble* xyz_list,GaussTria* gauss){ /*The Jacobian determinant is constant over the element, discard the gaussian points. * J is assumed to have been allocated of size NDOF2xNDOF2.*/ - double J[2][2]; + IssmDouble J[2][2]; /*Get Jacobian*/ GetJacobian(&J[0][0],xyz_list,gauss); @@ -357,11 +357,11 @@ } /*}}}*/ /*FUNCTION TriaRef::GetJacobianDeterminant3d {{{*/ -void TriaRef::GetJacobianDeterminant3d(double* Jdet, double* xyz_list,GaussTria* gauss){ +void TriaRef::GetJacobianDeterminant3d(IssmDouble* Jdet, IssmDouble* xyz_list,GaussTria* gauss){ /*The Jacobian determinant is constant over the element, discard the gaussian points. * J is assumed to have been allocated of size NDOF2xNDOF2.*/ - double x1,x2,x3,y1,y2,y3,z1,z2,z3; + IssmDouble x1,x2,x3,y1,y2,y3,z1,z2,z3; x1=*(xyz_list+3*0+0); y1=*(xyz_list+3*0+1); @@ -379,10 +379,10 @@ } /*}}}*/ /*FUNCTION TriaRef::GetJacobianInvert{{{*/ -void TriaRef::GetJacobianInvert(double* Jinv, double* xyz_list,GaussTria* gauss){ +void TriaRef::GetJacobianInvert(IssmDouble* Jinv, IssmDouble* xyz_list,GaussTria* gauss){ /*Jacobian*/ - double J[2][2]; + IssmDouble J[2][2]; /*Call Jacobian routine to get the jacobian:*/ GetJacobian(&J[0][0], xyz_list, gauss); @@ -393,7 +393,7 @@ } /*}}}*/ /*FUNCTION TriaRef::GetNodalFunctions{{{*/ -void TriaRef::GetNodalFunctions(double* basis,GaussTria* gauss){ +void TriaRef::GetNodalFunctions(IssmDouble* basis,GaussTria* gauss){ /*This routine returns the values of the nodal functions at the gaussian point.*/ basis[0]=gauss->coord1; @@ -403,10 +403,10 @@ } /*}}}*/ /*FUNCTION TriaRef::GetSegmentNodalFunctions{{{*/ -void TriaRef::GetSegmentNodalFunctions(double* basis,GaussTria* gauss,int index1,int index2){ +void TriaRef::GetSegmentNodalFunctions(IssmDouble* basis,GaussTria* gauss,int index1,int index2){ /*This routine returns the values of the nodal functions at the gaussian point.*/ - double BasisFunctions[3]; + IssmDouble BasisFunctions[3]; GetNodalFunctions(&BasisFunctions[0],gauss); @@ -417,13 +417,13 @@ } /*}}}*/ /*FUNCTION TriaRef::GetNodalFunctionsDerivatives{{{*/ -void TriaRef::GetNodalFunctionsDerivatives(double* dbasis,double* xyz_list, GaussTria* gauss){ +void TriaRef::GetNodalFunctionsDerivatives(IssmDouble* dbasis,IssmDouble* xyz_list, GaussTria* gauss){ /*This routine returns the values of the nodal functions derivatives (with respect to the * actual coordinate system): */ int i; - double dbasis_ref[NDOF2][NUMNODES]; - double Jinv[NDOF2][NDOF2]; + IssmDouble dbasis_ref[NDOF2][NUMNODES]; + IssmDouble Jinv[NDOF2][NDOF2]; /*Get derivative values with respect to parametric coordinate system: */ GetNodalFunctionsDerivativesReference(&dbasis_ref[0][0], gauss); @@ -444,7 +444,7 @@ } /*}}}*/ /*FUNCTION TriaRef::GetNodalFunctionsDerivativesReference{{{*/ -void TriaRef::GetNodalFunctionsDerivativesReference(double* dl1dl3,GaussTria* gauss){ +void TriaRef::GetNodalFunctionsDerivativesReference(IssmDouble* dl1dl3,GaussTria* gauss){ /*This routine returns the values of the nodal functions derivatives (with respect to the * natural coordinate system) at the gaussian point. */ @@ -463,7 +463,7 @@ } /*}}}*/ /*FUNCTION TriaRef::GetInputDerivativeValue{{{*/ -void TriaRef::GetInputDerivativeValue(double* p, double* plist,double* xyz_list, GaussTria* gauss){ +void TriaRef::GetInputDerivativeValue(IssmDouble* p, IssmDouble* plist,IssmDouble* xyz_list, GaussTria* gauss){ /*From node values of parameter p (plist[0],plist[1],plist[2]), return parameter derivative value at gaussian * point specified by gauss_basis: @@ -474,7 +474,7 @@ */ /*Nodal Derivatives*/ - double dbasis[2][3]; //nodal derivative functions in actual coordinate system. + IssmDouble dbasis[2][3]; //nodal derivative functions in actual coordinate system. /*Get dh1dh2dh3 in actual coordinate system: */ GetNodalFunctionsDerivatives(&dbasis[0][0],xyz_list, gauss); @@ -486,13 +486,13 @@ } /*}}}*/ /*FUNCTION TriaRef::GetInputValue{{{*/ -void TriaRef::GetInputValue(double* p, double* plist, GaussTria* gauss){ +void TriaRef::GetInputValue(IssmDouble* p, IssmDouble* plist, GaussTria* gauss){ /*From node values of parameter p (plist[0],plist[1],plist[2]), return parameter value at gaussian * point specifie by gauss: */ /*nodal functions annd output: */ - double basis[3]; + IssmDouble basis[3]; /*Get nodal functions*/ GetNodalFunctions(basis, gauss);