Index: /issm/trunk-jpl/src/c/objects/Elements/Element.h =================================================================== --- /issm/trunk-jpl/src/c/objects/Elements/Element.h (revision 12470) +++ /issm/trunk-jpl/src/c/objects/Elements/Element.h (revision 12471) @@ -38,11 +38,11 @@ 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; @@ -50,64 +50,64 @@ 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 }; Index: /issm/trunk-jpl/src/c/objects/Elements/Penta.cpp =================================================================== --- /issm/trunk-jpl/src/c/objects/Elements/Penta.cpp (revision 12470) +++ /issm/trunk-jpl/src/c/objects/Elements/Penta.cpp (revision 12471) @@ -145,15 +145,15 @@ /*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++){ @@ -181,6 +181,6 @@ /*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; @@ -233,16 +233,16 @@ 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; @@ -326,13 +326,13 @@ 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; @@ -764,9 +764,9 @@ /*}}}*/ /*FUNCTION Penta::GetElementSizes{{{*/ -void Penta::GetElementSizes(double* hx,double* hy,double* hz){ - - double xyz_list[NUMVERTICES][3]; - double xmin,ymin,zmin; - double xmax,ymax,zmax; +void Penta::GetElementSizes(IssmDouble* hx,IssmDouble* hy,IssmDouble* hz){ + + IssmDouble xyz_list[NUMVERTICES][3]; + IssmDouble xmin,ymin,zmin; + IssmDouble xmax,ymax,zmax; /*Get xyz list: */ @@ -820,9 +820,9 @@ } /*}}}*/ -/*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; @@ -845,9 +845,9 @@ } /*}}}*/ -/*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; @@ -874,6 +874,6 @@ } /*}}}*/ -/*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); @@ -888,10 +888,10 @@ /*}}}*/ /*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 */ @@ -978,11 +978,11 @@ /*}}}*/ /*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); @@ -996,5 +996,5 @@ /*}}}*/ /*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): * @@ -1008,6 +1008,6 @@ int i; - double epsilonvx[5]; - double epsilonvy[5]; + IssmDouble epsilonvx[5]; + IssmDouble epsilonvy[5]; /*Check that both inputs have been found*/ @@ -1025,5 +1025,5 @@ /*}}}*/ /*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): * @@ -1032,7 +1032,7 @@ 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*/ @@ -1100,10 +1100,10 @@ /*}}}*/ /*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); @@ -1127,5 +1127,5 @@ /*}}}*/ /*FUNCTION Penta::InputArtificialNoise{{{*/ -void Penta::InputArtificialNoise(int enum_type,double min,double max){ +void Penta::InputArtificialNoise(int enum_type,IssmDouble min,IssmDouble max){ Input* input=NULL; @@ -1140,5 +1140,5 @@ /*}}}*/ /*FUNCTION Penta::InputConvergence{{{*/ -bool Penta::InputConvergence(double* eps, int* enums,int num_enums,int* criterionenums,double* criterionvalues,int num_criterionenums){ +bool Penta::InputConvergence(IssmDouble* eps, int* enums,int num_enums,int* criterionenums,IssmDouble* criterionvalues,int num_criterionenums){ int i; @@ -1171,6 +1171,6 @@ } /*}}}*/ -/*FUNCTION Penta::InputCreate(double scalar,int enum,int code);{{{*/ -void Penta::InputCreate(double scalar,int name,int code){ +/*FUNCTION Penta::InputCreate(IssmDouble scalar,int enum,int code);{{{*/ +void Penta::InputCreate(IssmDouble scalar,int name,int code){ /*Check that name is an element input*/ @@ -1183,6 +1183,6 @@ this->inputs->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); @@ -1190,6 +1190,6 @@ } /*}}}*/ -/*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*/ @@ -1197,11 +1197,11 @@ 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: */ @@ -1222,5 +1222,5 @@ /*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: */ @@ -1237,5 +1237,5 @@ 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]; } @@ -1244,5 +1244,5 @@ /*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); @@ -1265,6 +1265,6 @@ 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); @@ -1284,7 +1284,7 @@ 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; @@ -1458,5 +1458,5 @@ /*}}}*/ /*FUNCTION Penta::InputScale{{{*/ -void Penta::InputScale(int enum_type,double scale_factor){ +void Penta::InputScale(int enum_type,IssmDouble scale_factor){ Input* input=NULL; @@ -1471,5 +1471,5 @@ /*}}}*/ /*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; @@ -1505,6 +1505,6 @@ } /*}}}*/ -/*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; @@ -1529,9 +1529,9 @@ IssmInt i,j; int penta_vertex_ids[6]; - double nodeinputs[6]; - double cmmininputs[6]; - double cmmaxinputs[6]; - - double yts; + IssmDouble nodeinputs[6]; + IssmDouble cmmininputs[6]; + IssmDouble cmmaxinputs[6]; + + IssmDouble yts; bool control_analysis; int num_control_type; @@ -1647,5 +1647,5 @@ /*}}}*/ /*FUNCTION Penta::InputUpdateFromSolution {{{*/ -void Penta::InputUpdateFromSolution(double* solution){ +void Penta::InputUpdateFromSolution(IssmDouble* solution){ int analysis_type; @@ -1716,5 +1716,5 @@ /*}}}*/ /*FUNCTION Penta::InputUpdateFromSolutionPrognostic{{{*/ -void Penta::InputUpdateFromSolutionPrognostic(double* solution){ +void Penta::InputUpdateFromSolutionPrognostic(IssmDouble* solution){ const int numdof = NDOF1*NUMVERTICES; @@ -1723,11 +1723,11 @@ 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; @@ -1799,9 +1799,9 @@ /*}}}*/ /*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; @@ -1823,10 +1823,10 @@ /*}}}*/ /*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; @@ -1862,6 +1862,6 @@ } /*}}}*/ -/*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*/ @@ -1875,5 +1875,5 @@ /*New PentaVertexInpu*/ - double values[6]; + IssmDouble values[6]; /*Get values on the 6 vertices*/ @@ -1978,5 +1978,5 @@ /*}}}*/ /*FUNCTION Penta::IsNodeOnShelfFromFlags {{{*/ -bool Penta::IsNodeOnShelfFromFlags(double* flags){ +bool Penta::IsNodeOnShelfFromFlags(IssmDouble* flags){ int i; @@ -2009,5 +2009,5 @@ /*}}}*/ /*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*/ @@ -2016,5 +2016,5 @@ int *resultsenums = NULL; int *resultssizes = NULL; - double *resultstimes = NULL; + IssmDouble *resultstimes = NULL; int *resultssteps = NULL; @@ -2033,5 +2033,5 @@ resultsenums=xNew(numberofresults); resultssizes=xNew(numberofresults); - resultstimes=xNew(numberofresults); + resultstimes=xNew(numberofresults); resultssteps=xNew(numberofresults); @@ -2060,12 +2060,12 @@ }/*}}}*/ /*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; @@ -2143,11 +2143,11 @@ /*}}}*/ /*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++){ @@ -2171,9 +2171,9 @@ /*}}}*/ /*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; @@ -2253,58 +2253,58 @@ /*}}}*/ /*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 - - 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 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 + + 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}; - - double tstar; // monthly mean surface temp - double Tsum[NUMVERTICES]= {0}; // average summer (JJA) temperature - double Tsurf[NUMVERTICES] = {0}; // average annual temperature + 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}; + + 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). - - 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 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). + + 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: */ @@ -2481,7 +2481,7 @@ 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; @@ -2518,13 +2518,13 @@ /*}}}*/ /*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 */ @@ -2560,13 +2560,13 @@ /*}}}*/ /*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 */ @@ -2595,5 +2595,5 @@ /*}}}*/ /*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)){ @@ -2649,7 +2649,7 @@ 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*/ @@ -2739,8 +2739,8 @@ /*}}}*/ /*FUNCTION Penta::SurfaceArea {{{*/ -double Penta::SurfaceArea(void){ +IssmDouble Penta::SurfaceArea(void){ int approximation; - double S; + IssmDouble S; Tria* tria=NULL; @@ -2776,10 +2776,10 @@ /*}}}*/ /*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++){ @@ -2800,11 +2800,11 @@ /*}}}*/ /*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:*/ @@ -2851,10 +2851,10 @@ 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: */ @@ -2998,5 +2998,5 @@ /*}}}*/ /*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; @@ -3025,10 +3025,10 @@ /*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; @@ -3066,5 +3066,5 @@ /*}}}*/ /*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"); } @@ -3073,9 +3073,9 @@ #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; @@ -3096,8 +3096,8 @@ /*}}}*/ /*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: */ @@ -3109,8 +3109,8 @@ /*}}}*/ /*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: */ @@ -3122,8 +3122,8 @@ /*}}}*/ /*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: */ @@ -3135,8 +3135,8 @@ /*}}}*/ /*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: */ @@ -3148,7 +3148,7 @@ /*}}}*/ /*FUNCTION Penta::MassFlux {{{*/ -double Penta::MassFlux( double* segment,bool process_units){ - - double mass_flux=0; +IssmDouble Penta::MassFlux( IssmDouble* segment,bool process_units){ + + IssmDouble mass_flux=0; if(!IsOnBed()) return mass_flux; @@ -3172,8 +3172,8 @@ /*}}}*/ /*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: */ @@ -3185,8 +3185,8 @@ /*}}}*/ /*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: */ @@ -3198,8 +3198,8 @@ /*}}}*/ /*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: */ @@ -3211,8 +3211,8 @@ /*}}}*/ /*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: */ @@ -3225,8 +3225,8 @@ /*}}}*/ /*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: */ @@ -3238,8 +3238,8 @@ /*}}}*/ /*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: */ @@ -3251,8 +3251,8 @@ /*}}}*/ /*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: */ @@ -3264,5 +3264,5 @@ /*}}}*/ /*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){ @@ -3273,5 +3273,5 @@ /*Get input:*/ - double vel; + IssmDouble vel; Input* vel_input; @@ -3316,22 +3316,22 @@ 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; @@ -3472,11 +3472,11 @@ /*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; @@ -3553,19 +3553,19 @@ 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; @@ -3705,11 +3705,11 @@ /*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; @@ -3776,17 +3776,17 @@ 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; @@ -3871,11 +3871,11 @@ /*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; @@ -3929,12 +3929,12 @@ 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; @@ -4039,16 +4039,16 @@ 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; @@ -4125,11 +4125,11 @@ /*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; @@ -4183,11 +4183,11 @@ 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; @@ -4248,6 +4248,6 @@ int i; int* doflist=NULL; - double values[numdof]; - double temp; + IssmDouble values[numdof]; + IssmDouble temp; GaussPenta *gauss=NULL; @@ -4279,6 +4279,6 @@ int i; int* doflist=NULL; - double values[numdof]; - double enthalpy; + IssmDouble values[numdof]; + IssmDouble enthalpy; GaussPenta *gauss=NULL; @@ -4304,5 +4304,5 @@ /*}}}*/ /*FUNCTION Penta::InputUpdateFromSolutionThermal {{{*/ -void Penta::InputUpdateFromSolutionThermal(double* solution){ +void Penta::InputUpdateFromSolutionThermal(IssmDouble* solution){ const int numdof=NDOF1*NUMVERTICES; @@ -4310,10 +4310,10 @@ 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: */ @@ -4375,5 +4375,5 @@ /*}}}*/ /*FUNCTION Penta::InputUpdateFromSolutionEnthalpy {{{*/ -void Penta::InputUpdateFromSolutionEnthalpy(double* solution){ +void Penta::InputUpdateFromSolutionEnthalpy(IssmDouble* solution){ const int numdof=NDOF1*NUMVERTICES; @@ -4381,11 +4381,11 @@ 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; @@ -4476,5 +4476,5 @@ }/*}}}*/ /*FUNCTION Penta::ControlInputScaleGradient{{{*/ -void Penta::ControlInputScaleGradient(int enum_type,double scale){ +void Penta::ControlInputScaleGradient(int enum_type,IssmDouble scale){ Input* input=NULL; @@ -4492,8 +4492,8 @@ }/*}}}*/ /*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; @@ -4566,13 +4566,13 @@ 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; @@ -4634,14 +4634,14 @@ 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; @@ -4889,12 +4889,12 @@ 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; @@ -4960,14 +4960,14 @@ 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; @@ -5099,5 +5099,5 @@ } /*}}}*/ /*FUNCTION Penta::InputControlUpdate{{{*/ -void Penta::InputControlUpdate(double scalar,bool save_parameter){ +void Penta::InputControlUpdate(IssmDouble scalar,bool save_parameter){ /*Intermediary*/ @@ -5137,14 +5137,14 @@ /*}}}*/ /*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; @@ -5180,12 +5180,12 @@ /*}}}*/ /*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; @@ -5215,8 +5215,8 @@ /*}}}*/ /*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; @@ -5252,8 +5252,8 @@ /*}}}*/ /*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; @@ -5289,8 +5289,8 @@ /*}}}*/ /*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; @@ -5326,7 +5326,7 @@ /*}}}*/ /*FUNCTION Penta::SurfaceLogVxVyMisfit {{{*/ -double Penta::SurfaceLogVxVyMisfit(bool process_units,int weight_index){ - - double J; +IssmDouble Penta::SurfaceLogVxVyMisfit(bool process_units,int weight_index){ + + IssmDouble J; Tria* tria=NULL; @@ -5365,8 +5365,8 @@ /*}}}*/ /*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; @@ -5402,5 +5402,5 @@ /*}}}*/ /*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"); @@ -5408,8 +5408,8 @@ /*}}}*/ /*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; @@ -5428,7 +5428,7 @@ /*}}}*/ /*FUNCTION Penta::DragCoefficientAbsGradient{{{*/ -double Penta::DragCoefficientAbsGradient(bool process_units,int weight_index){ - - double J; +IssmDouble Penta::DragCoefficientAbsGradient(bool process_units,int weight_index){ + + IssmDouble J; Tria* tria=NULL; @@ -5443,7 +5443,7 @@ /*}}}*/ /*FUNCTION Penta::RheologyBbarAbsGradient{{{*/ -double Penta::RheologyBbarAbsGradient(bool process_units,int weight_index){ - - double J; +IssmDouble Penta::RheologyBbarAbsGradient(bool process_units,int weight_index){ + + IssmDouble J; Tria* tria=NULL; @@ -5481,7 +5481,7 @@ /*}}}*/ /*FUNCTION Penta::SetControlInputsFromVector{{{*/ -void Penta::SetControlInputsFromVector(double* vector,int control_enum,int control_index){ - - double values[NUMVERTICES]; +void Penta::SetControlInputsFromVector(IssmDouble* vector,int control_enum,int control_index){ + + IssmDouble values[NUMVERTICES]; int doflist1[NUMVERTICES]; Input *input = NULL; @@ -5518,6 +5518,6 @@ #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; @@ -5531,5 +5531,5 @@ /*New PentaP1Input*/ - double values[6]; + IssmDouble values[6]; /*Get values on the 6 vertices*/ @@ -5542,9 +5542,9 @@ 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: */ @@ -5560,5 +5560,5 @@ 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(); @@ -5630,13 +5630,13 @@ } /*}}}*/ -/*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*/ @@ -5656,9 +5656,9 @@ 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); @@ -5742,14 +5742,14 @@ /*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; @@ -5834,15 +5834,15 @@ /*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; @@ -5892,5 +5892,5 @@ if (slope_magnitude>MAXSLOPE){ - alpha2=pow((double)10,MOUNTAINKEXPONENT); + alpha2=pow((IssmDouble)10,MOUNTAINKEXPONENT); } @@ -5947,19 +5947,19 @@ /*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; @@ -6054,18 +6054,18 @@ 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; @@ -6242,5 +6242,5 @@ int connectivity[2]; int i,i0,i1,j0,j1; - double one0,one1; + IssmDouble one0,one1; /*Initialize Element matrix*/ @@ -6259,6 +6259,6 @@ 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*/ @@ -6341,14 +6341,14 @@ /*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; @@ -6494,12 +6494,12 @@ 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; @@ -6559,13 +6559,13 @@ 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; @@ -6604,5 +6604,5 @@ //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); } @@ -6660,12 +6660,12 @@ /*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; @@ -6726,12 +6726,12 @@ 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; @@ -6811,9 +6811,9 @@ /*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; @@ -6857,9 +6857,9 @@ /*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; @@ -6917,11 +6917,11 @@ 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; @@ -6979,13 +6979,13 @@ 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; @@ -7068,11 +7068,11 @@ 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; @@ -7130,13 +7130,13 @@ 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; @@ -7280,12 +7280,12 @@ 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; @@ -7335,8 +7335,8 @@ 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]; } } @@ -7345,10 +7345,10 @@ //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]; } } @@ -7380,9 +7380,9 @@ /*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; @@ -7443,17 +7443,17 @@ 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; @@ -7520,12 +7520,12 @@ 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; @@ -7605,9 +7605,9 @@ 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; @@ -7662,10 +7662,10 @@ 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; @@ -7768,13 +7768,13 @@ /*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; @@ -7833,14 +7833,14 @@ /*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; @@ -7909,6 +7909,6 @@ int approximation; int* doflist=NULL; - double vx,vy; - double values[numdof]; + IssmDouble vx,vy; + IssmDouble values[numdof]; GaussPenta* gauss; @@ -7950,6 +7950,6 @@ int i; int* doflist=NULL; - double vx,vy; - double values[numdof]; + IssmDouble vx,vy; + IssmDouble values[numdof]; GaussPenta* gauss=NULL; @@ -7986,6 +7986,6 @@ int i; int* doflist=NULL; - double vz; - double values[numdof]; + IssmDouble vz; + IssmDouble values[numdof]; GaussPenta* gauss=NULL; @@ -8019,7 +8019,7 @@ 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; @@ -8058,5 +8058,5 @@ /*}}}*/ /*FUNCTION Penta::InputUpdateFromSolutionDiagnosticHoriz {{{*/ -void Penta::InputUpdateFromSolutionDiagnosticHoriz(double* solution){ +void Penta::InputUpdateFromSolutionDiagnosticHoriz(IssmDouble* solution){ int approximation; @@ -8094,18 +8094,18 @@ /*}}}*/ /*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; @@ -8175,5 +8175,5 @@ /*}}}*/ /*FUNCTION Penta::InputUpdateFromSolutionDiagnosticMacAyealPattyn {{{*/ -void Penta::InputUpdateFromSolutionDiagnosticMacAyealPattyn(double* solution){ +void Penta::InputUpdateFromSolutionDiagnosticMacAyealPattyn(IssmDouble* solution){ const int numdof=NDOF2*NUMVERTICES; @@ -8181,14 +8181,14 @@ 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; @@ -8259,5 +8259,5 @@ /*}}}*/ /*FUNCTION Penta::InputUpdateFromSolutionDiagnosticMacAyealStokes {{{*/ -void Penta::InputUpdateFromSolutionDiagnosticMacAyealStokes(double* solution){ +void Penta::InputUpdateFromSolutionDiagnosticMacAyealStokes(IssmDouble* solution){ const int numdofm=NDOF2*NUMVERTICES; @@ -8266,15 +8266,15 @@ 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; @@ -8359,18 +8359,18 @@ /*}}}*/ /*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; @@ -8433,5 +8433,5 @@ /*}}}*/ /*FUNCTION Penta::InputUpdateFromSolutionDiagnosticPattynStokes {{{*/ -void Penta::InputUpdateFromSolutionDiagnosticPattynStokes(double* solution){ +void Penta::InputUpdateFromSolutionDiagnosticPattynStokes(IssmDouble* solution){ const int numdofp=NDOF2*NUMVERTICES; @@ -8439,15 +8439,15 @@ 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; @@ -8527,18 +8527,18 @@ /*}}}*/ /*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; @@ -8590,5 +8590,5 @@ /*}}}*/ /*FUNCTION Penta::InputUpdateFromSolutionDiagnosticVert {{{*/ -void Penta::InputUpdateFromSolutionDiagnosticVert(double* solution){ +void Penta::InputUpdateFromSolutionDiagnosticVert(IssmDouble* solution){ const int numdof=NDOF1*NUMVERTICES; @@ -8596,16 +8596,16 @@ 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; @@ -8694,16 +8694,16 @@ /*}}}*/ /*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; Index: /issm/trunk-jpl/src/c/objects/Elements/Penta.h =================================================================== --- /issm/trunk-jpl/src/c/objects/Elements/Penta.h (revision 12470) +++ /issm/trunk-jpl/src/c/objects/Elements/Penta.h (revision 12471) @@ -59,20 +59,20 @@ /*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); @@ -87,55 +87,55 @@ 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); @@ -147,21 +147,21 @@ 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); @@ -173,20 +173,20 @@ 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); @@ -194,12 +194,12 @@ 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_ @@ -236,13 +236,13 @@ 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); @@ -278,6 +278,6 @@ ElementVector* CreatePVectorAdjointPattyn(void); ElementVector* CreatePVectorAdjointStokes(void); - void InputUpdateFromSolutionAdjointHoriz( double* solutiong); - void InputUpdateFromSolutionAdjointStokes( double* solutiong); + void InputUpdateFromSolutionAdjointHoriz( IssmDouble* solutiong); + void InputUpdateFromSolutionAdjointStokes( IssmDouble* solutiong); #endif @@ -303,6 +303,6 @@ 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_ Index: /issm/trunk-jpl/src/c/objects/Elements/PentaRef.cpp =================================================================== --- /issm/trunk-jpl/src/c/objects/Elements/PentaRef.cpp (revision 12470) +++ /issm/trunk-jpl/src/c/objects/Elements/PentaRef.cpp (revision 12471) @@ -58,5 +58,5 @@ /*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 @@ -70,5 +70,5 @@ */ - double dbasis[3][NUMNODESP1]; + IssmDouble dbasis[3][NUMNODESP1]; /*Get dbasis in actual coordinate system: */ @@ -89,5 +89,5 @@ /*}}}*/ /*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 @@ -103,6 +103,6 @@ int i; - double dh1dh7[3][NUMNODESMINI]; - double l1l6[NUMNODESP1]; + IssmDouble dh1dh7[3][NUMNODESMINI]; + IssmDouble l1l6[NUMNODESP1]; /*Get dh1dh6 in actual coordinate system: */ @@ -135,5 +135,5 @@ /*}}}*/ /*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 @@ -149,5 +149,5 @@ */ - double dbasis[3][NUMNODESP1]; + IssmDouble dbasis[3][NUMNODESP1]; /*Get dbasis in actual coordinate system: */ @@ -175,5 +175,5 @@ /*}}}*/ /*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 @@ -188,5 +188,5 @@ * 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: */ @@ -213,5 +213,5 @@ /*}}}*/ /*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 @@ -226,5 +226,5 @@ int i; - double dh1dh7[3][NUMNODESMINI]; + IssmDouble dh1dh7[3][NUMNODESMINI]; /*Get dh1dh6 in actual coordinate system: */ @@ -253,5 +253,5 @@ /*}}}*/ /*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. @@ -271,6 +271,6 @@ int i; - double dh1dh7[3][NUMNODESMINI]; - double l1l6[NUMNODESP1]; + IssmDouble dh1dh7[3][NUMNODESMINI]; + IssmDouble l1l6[NUMNODESP1]; /*Get dh1dh7 in actual coordinate system: */ @@ -320,5 +320,5 @@ /*}}}*/ /*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 @@ -338,6 +338,6 @@ int i; - double dh1dh7[3][NUMNODESMINI]; - double l1l6[NUMNODESP1]; + IssmDouble dh1dh7[3][NUMNODESMINI]; + IssmDouble l1l6[NUMNODESP1]; /*Get dh1dh7 in actual coordinate system: */ @@ -387,5 +387,5 @@ /*}}}*/ /*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 @@ -400,5 +400,5 @@ /*Same thing in the actual coordinate system: */ - double l1l6[6]; + IssmDouble l1l6[6]; /*Get dh1dh2dh3 in actual coordinates system : */ @@ -414,5 +414,5 @@ /*}}}*/ /*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 @@ -427,5 +427,5 @@ /*Same thing in the actual coordinate system: */ - double dh1dh6[3][NUMNODESP1]; + IssmDouble dh1dh6[3][NUMNODESP1]; /*Get dh1dh2dh3 in actual coordinates system : */ @@ -441,10 +441,10 @@ /*}}}*/ /*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: */ @@ -459,5 +459,5 @@ /*}}}*/ /*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 @@ -472,5 +472,5 @@ /*Same thing in the actual coordinate system: */ - double dh1dh6[3][NUMNODESP1]; + IssmDouble dh1dh6[3][NUMNODESP1]; /*Get dh1dh2dh3 in actual coordinates system : */ @@ -486,5 +486,5 @@ /*}}}*/ /*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*/ @@ -494,5 +494,5 @@ /*}}}*/ /*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 @@ -509,5 +509,5 @@ int i; - double l1l6[6]; + IssmDouble l1l6[6]; /*Get l1l6 in actual coordinate system: */ @@ -531,5 +531,5 @@ /*}}}*/ /*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. @@ -544,5 +544,5 @@ const int num_dof=4; - double l1l2l3[NUMNODESP1_2d]; + IssmDouble l1l2l3[NUMNODESP1_2d]; /*Get l1l2l3 in actual coordinate system: */ @@ -566,5 +566,5 @@ /*}}}*/ /*FUNCTION PentaRef::GetLprimeStokes {{{*/ -void PentaRef::GetLprimeStokes(double* LprimeStokes, double* xyz_list, GaussPenta* gauss){ +void PentaRef::GetLprimeStokes(IssmDouble* LprimeStokes, IssmDouble* xyz_list, GaussPenta* gauss){ /* @@ -606,6 +606,6 @@ 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: */ @@ -678,5 +678,5 @@ /*}}}*/ /*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. @@ -697,5 +697,5 @@ int num_dof=2; - double l1l2l3[NUMNODESP1_2d]; + IssmDouble l1l2l3[NUMNODESP1_2d]; @@ -728,5 +728,5 @@ /*}}}*/ /*FUNCTION PentaRef::GetLprimeMacAyealStokes {{{*/ -void PentaRef::GetLprimeMacAyealStokes(double* LprimeStokes, double* xyz_list, GaussPenta* gauss){ +void PentaRef::GetLprimeMacAyealStokes(IssmDouble* LprimeStokes, IssmDouble* xyz_list, GaussPenta* gauss){ /* @@ -747,6 +747,6 @@ 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: */ @@ -795,5 +795,5 @@ /*}}}*/ /*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. @@ -810,5 +810,5 @@ int num_dof=4; - double l1l2l3[NUMNODESP1_2d]; + IssmDouble l1l2l3[NUMNODESP1_2d]; @@ -841,5 +841,5 @@ /*}}}*/ /*FUNCTION PentaRef::GetLprimeStokesMacAyeal {{{*/ -void PentaRef::GetLprimeStokesMacAyeal(double* LprimeStokes, double* xyz_list, GaussPenta* gauss){ +void PentaRef::GetLprimeStokesMacAyeal(IssmDouble* LprimeStokes, IssmDouble* xyz_list, GaussPenta* gauss){ /* @@ -856,6 +856,6 @@ 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: */ @@ -880,5 +880,5 @@ /*}}}*/ /*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; @@ -887,10 +887,10 @@ * J is assumed to have been allocated of size NDOF2xNDOF2.*/ - double A1,A2,A3; //area coordinates - double 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 A1,A2,A3; //area coordinates + IssmDouble xi,eta,zi; //parametric coordinates + + 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: */ @@ -939,8 +939,8 @@ /*}}}*/ /*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*/ @@ -954,9 +954,9 @@ /*}}}*/ /*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); @@ -976,9 +976,9 @@ /*}}}*/ /*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); @@ -995,8 +995,8 @@ /*}}}*/ /*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:*/ @@ -1008,5 +1008,5 @@ /*}}}*/ /*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.*/ @@ -1022,5 +1022,5 @@ /*}}}*/ /*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 @@ -1028,6 +1028,6 @@ 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: */ @@ -1053,11 +1053,11 @@ /*}}}*/ /*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: */ @@ -1099,5 +1099,5 @@ /*}}}*/ /*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.*/ @@ -1112,10 +1112,10 @@ /*}}}*/ /*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: */ @@ -1141,10 +1141,10 @@ /*}}}*/ /*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); @@ -1185,8 +1185,8 @@ /*}}}*/ /*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); @@ -1205,8 +1205,8 @@ /*}}}*/ /*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]; @@ -1231,9 +1231,9 @@ /*}}}*/ /*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: */ @@ -1246,5 +1246,5 @@ /*}}}*/ /*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; @@ -1254,5 +1254,5 @@ * 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: */ Index: /issm/trunk-jpl/src/c/objects/Elements/PentaRef.h =================================================================== --- /issm/trunk-jpl/src/c/objects/Elements/PentaRef.h (revision 12470) +++ /issm/trunk-jpl/src/c/objects/Elements/PentaRef.h (revision 12471) @@ -23,40 +23,40 @@ /*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");}; }; Index: /issm/trunk-jpl/src/c/objects/Elements/Tria.cpp =================================================================== --- /issm/trunk-jpl/src/c/objects/Elements/Tria.cpp (revision 12470) +++ /issm/trunk-jpl/src/c/objects/Elements/Tria.cpp (revision 12471) @@ -125,5 +125,5 @@ /*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; @@ -132,7 +132,7 @@ int offsetsid[NUMVERTICES]; int offsetdof[NUMVERTICES]; - double area; - double mean; - double values[3]; + IssmDouble area; + IssmDouble mean; + IssmDouble values[3]; /*First, get the area: */ @@ -232,8 +232,8 @@ /*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; @@ -291,16 +291,16 @@ 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; @@ -410,12 +410,12 @@ /*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; @@ -485,7 +485,7 @@ /* 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; @@ -599,8 +599,8 @@ /*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; @@ -650,8 +650,8 @@ /*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; @@ -696,8 +696,8 @@ 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; @@ -787,13 +787,13 @@ 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; @@ -948,9 +948,9 @@ /*}}}*/ /*FUNCTION Tria::GetArea {{{*/ -double Tria::GetArea(void){ - - double area=0; - double xyz_list[NUMVERTICES][3]; - double x1,y1,x2,y2,x3,y3; +IssmDouble Tria::GetArea(void){ + + IssmDouble area=0; + IssmDouble xyz_list[NUMVERTICES][3]; + IssmDouble x1,y1,x2,y2,x3,y3; /*Get xyz list: */ @@ -1022,9 +1022,9 @@ } /*}}}*/ -/*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; @@ -1047,8 +1047,8 @@ } /*}}}*/ -/*FUNCTION Tria::GetInputListOnVertices(double* pvalue,int enumtype,double defaultvalue) {{{*/ -void Tria::GetInputListOnVertices(double* pvalue,int enumtype,double defaultvalue){ - - double value[NUMVERTICES]; +/*FUNCTION Tria::GetInputListOnVertices(IssmDouble* pvalue,int enumtype,IssmDouble defaultvalue) {{{*/ +void Tria::GetInputListOnVertices(IssmDouble* pvalue,int enumtype,IssmDouble defaultvalue){ + + IssmDouble value[NUMVERTICES]; GaussTria *gauss = NULL; Input *input = inputs->GetInput(enumtype); @@ -1073,8 +1073,8 @@ } /*}}}*/ -/*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){ - - double value[NUMVERTICES]; +/*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){ + + IssmDouble value[NUMVERTICES]; GaussTria *gauss = NULL; Input *input = inputs->GetInput(enumtype); @@ -1099,6 +1099,6 @@ } /*}}}*/ -/*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); @@ -1150,12 +1150,12 @@ } /*}}}*/ -/*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*/ @@ -1226,5 +1226,5 @@ /*}}}*/ /*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; @@ -1239,5 +1239,5 @@ /*}}}*/ /*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; @@ -1306,5 +1306,5 @@ /*}}}*/ /*FUNCTION Tria::InputScale{{{*/ -void Tria::InputScale(int enum_type,double scale_factor){ +void Tria::InputScale(int enum_type,IssmDouble scale_factor){ Input* input=NULL; @@ -1319,5 +1319,5 @@ /*}}}*/ /*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; @@ -1350,6 +1350,6 @@ } /*}}}*/ -/*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; @@ -1374,10 +1374,10 @@ 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; @@ -1455,5 +1455,5 @@ /*}}}*/ /*FUNCTION Tria::InputUpdateFromSolution {{{*/ -void Tria::InputUpdateFromSolution(double* solution){ +void Tria::InputUpdateFromSolution(IssmDouble* solution){ /*retrive parameters: */ @@ -1510,10 +1510,10 @@ /*}}}*/ /*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: */ @@ -1534,5 +1534,5 @@ /*}}}*/ /*FUNCTION Tria::InputUpdateFromSolutionPrognostic{{{*/ -void Tria::InputUpdateFromSolutionPrognostic(double* solution){ +void Tria::InputUpdateFromSolutionPrognostic(IssmDouble* solution){ /*Intermediaries*/ @@ -1541,11 +1541,11 @@ 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: */ @@ -1600,6 +1600,6 @@ } /*}}}*/ -/*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*/ @@ -1611,5 +1611,5 @@ /*New TriaP1Input*/ - double values[3]; + IssmDouble values[3]; /*Get values on the 3 vertices*/ @@ -1642,6 +1642,6 @@ } /*}}}*/ -/*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*/ @@ -1654,5 +1654,5 @@ 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)); } @@ -1661,6 +1661,6 @@ } /*}}}*/ -/*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*/ @@ -1668,10 +1668,10 @@ 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; @@ -1693,5 +1693,5 @@ /*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: */ @@ -1708,5 +1708,5 @@ 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]; } @@ -1715,5 +1715,5 @@ /*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); @@ -1736,6 +1736,6 @@ 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); @@ -1817,5 +1817,5 @@ /*}}}*/ /*FUNCTION Tria::IsNodeOnShelfFromFlags {{{*/ -bool Tria::IsNodeOnShelfFromFlags(double* flags){ +bool Tria::IsNodeOnShelfFromFlags(IssmDouble* flags){ int i; @@ -1840,5 +1840,5 @@ /*}}}*/ /*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*/ @@ -1847,5 +1847,5 @@ int *resultsenums = NULL; int *resultssizes = NULL; - double *resultstimes = NULL; + IssmDouble *resultstimes = NULL; int *resultssteps = NULL; @@ -1864,5 +1864,5 @@ resultsenums=xNew(numberofresults); resultssizes=xNew(numberofresults); - resultstimes=xNew(numberofresults); + resultstimes=xNew(numberofresults); resultssteps=xNew(numberofresults); @@ -1891,12 +1891,12 @@ }/*}}}*/ /*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: */ @@ -1973,9 +1973,9 @@ /*}}}*/ /*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; @@ -2058,7 +2058,7 @@ 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; @@ -2084,56 +2084,56 @@ /*}}}*/ /*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 - - 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 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 + + 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}; - - double tstar; // monthly mean surface temp - double Tsum[NUMVERTICES]= {0}; // average summer (JJA) temperature - double Tsurf[NUMVERTICES] = {0}; // average annual temperature + 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}; + + 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). - - 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 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). + + 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: */ @@ -2146,5 +2146,5 @@ 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); @@ -2330,5 +2330,5 @@ /*}}}*/ /*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)){ @@ -2364,5 +2364,5 @@ /*}}}*/ /*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"); @@ -2386,11 +2386,11 @@ /*}}}*/ /*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: */ @@ -2408,5 +2408,5 @@ 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: */ @@ -2415,10 +2415,10 @@ /*}}}*/ /*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++){ @@ -2431,5 +2431,5 @@ 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; @@ -2439,14 +2439,14 @@ /*}}}*/ /*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:*/ @@ -2492,6 +2492,6 @@ int tria_vertex_ids[3]; int tria_type; - double nodeinputs[3]; - double yts; + IssmDouble nodeinputs[3]; + IssmDouble yts; int progstabilization,balancestabilization; bool dakota_analysis; @@ -2585,5 +2585,5 @@ /*}}}*/ /*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; @@ -2607,9 +2607,9 @@ #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; @@ -2633,15 +2633,15 @@ /*}}}*/ /*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; @@ -2705,8 +2705,8 @@ /*}}}*/ /*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: */ @@ -2718,8 +2718,8 @@ /*}}}*/ /*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: */ @@ -2731,8 +2731,8 @@ /*}}}*/ /*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: */ @@ -2744,8 +2744,8 @@ /*}}}*/ /*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: */ @@ -2757,8 +2757,8 @@ /*}}}*/ /*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: */ @@ -2770,8 +2770,8 @@ /*}}}*/ /*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: */ @@ -2784,8 +2784,8 @@ /*}}}*/ /*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: */ @@ -2797,8 +2797,8 @@ /*}}}*/ /*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: */ @@ -2810,8 +2810,8 @@ /*}}}*/ /*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: */ @@ -2823,8 +2823,8 @@ /*}}}*/ /*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: */ @@ -2836,8 +2836,8 @@ /*}}}*/ /*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: */ @@ -2849,5 +2849,5 @@ /*}}}*/ /*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){ @@ -2858,5 +2858,5 @@ /*Get input:*/ - double vel; + IssmDouble vel; Input* vel_input; @@ -2900,12 +2900,12 @@ /*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; @@ -2965,13 +2965,13 @@ 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; @@ -3002,5 +3002,5 @@ 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); @@ -3038,6 +3038,6 @@ 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; } @@ -3054,10 +3054,10 @@ /*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; @@ -3105,8 +3105,8 @@ /*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; @@ -3138,9 +3138,9 @@ 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]; - - 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; + 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])/(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; } @@ -3158,13 +3158,13 @@ /*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; @@ -3223,6 +3223,6 @@ int i; int* doflist=NULL; - double vx,vy; - double values[numdof]; + IssmDouble vx,vy; + IssmDouble values[numdof]; GaussTria* gauss=NULL; @@ -3261,6 +3261,6 @@ int i; - double vx,vy; - double values[numdof]; + IssmDouble vx,vy; + IssmDouble values[numdof]; int *doflist = NULL; GaussTria *gauss = NULL; @@ -3295,5 +3295,5 @@ /*}}}*/ /*FUNCTION Tria::InputUpdateFromSolutionDiagnosticHoriz {{{*/ -void Tria::InputUpdateFromSolutionDiagnosticHoriz(double* solution){ +void Tria::InputUpdateFromSolutionDiagnosticHoriz(IssmDouble* solution){ const int numdof=NDOF2*NUMVERTICES; @@ -3301,12 +3301,12 @@ 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: */ @@ -3358,5 +3358,5 @@ /*}}}*/ /*FUNCTION Tria::InputUpdateFromSolutionDiagnosticHutter {{{*/ -void Tria::InputUpdateFromSolutionDiagnosticHutter(double* solution){ +void Tria::InputUpdateFromSolutionDiagnosticHutter(IssmDouble* solution){ const int numdof=NDOF2*NUMVERTICES; @@ -3364,12 +3364,12 @@ 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: */ @@ -3420,5 +3420,5 @@ #ifdef _HAVE_CONTROL_ /*FUNCTION Tria::InputControlUpdate{{{*/ -void Tria::InputControlUpdate(double scalar,bool save_parameter){ +void Tria::InputControlUpdate(IssmDouble scalar,bool save_parameter){ /*Intermediary*/ @@ -3474,5 +3474,5 @@ }/*}}}*/ /*FUNCTION Tria::ControlInputScaleGradient{{{*/ -void Tria::ControlInputScaleGradient(int enum_type,double scale){ +void Tria::ControlInputScaleGradient(int enum_type,IssmDouble scale){ Input* input=NULL; @@ -3490,8 +3490,8 @@ }/*}}}*/ /*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; @@ -3577,9 +3577,9 @@ 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; @@ -3618,10 +3618,10 @@ 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; @@ -3676,12 +3676,12 @@ 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; @@ -3746,5 +3746,5 @@ // 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*/ @@ -3762,9 +3762,9 @@ 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; @@ -3806,6 +3806,6 @@ /*Intermediaries*/ int doflist1[NUMVERTICES]; - double lambda[NUMVERTICES]; - double gradient_g[NUMVERTICES]; + IssmDouble lambda[NUMVERTICES]; + IssmDouble gradient_g[NUMVERTICES]; /*Compute Gradient*/ @@ -3823,9 +3823,9 @@ 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; @@ -3866,9 +3866,9 @@ 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; @@ -3917,13 +3917,13 @@ /*}}}*/ /*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; @@ -3961,13 +3961,13 @@ /*}}}*/ /*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; @@ -4022,16 +4022,16 @@ /*}}}*/ /*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; @@ -4087,5 +4087,5 @@ /*}}}*/ /*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; @@ -4093,10 +4093,10 @@ 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; @@ -4153,13 +4153,13 @@ /*}}}*/ /*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; @@ -4214,15 +4214,15 @@ /*}}}*/ /*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; @@ -4278,13 +4278,13 @@ /*}}}*/ /*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; @@ -4322,14 +4322,14 @@ /*}}}*/ /*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: */ @@ -4374,12 +4374,12 @@ /*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; @@ -4442,13 +4442,13 @@ 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; @@ -4618,13 +4618,13 @@ 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; @@ -4789,13 +4789,13 @@ /*}}}*/ /*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; @@ -4863,13 +4863,13 @@ 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; @@ -4925,5 +4925,5 @@ /*}}}*/ /*FUNCTION Tria::InputUpdateFromSolutionAdjointHoriz {{{*/ -void Tria::InputUpdateFromSolutionAdjointHoriz(double* solution){ +void Tria::InputUpdateFromSolutionAdjointHoriz(IssmDouble* solution){ const int numdof=NDOF2*NUMVERTICES; @@ -4931,7 +4931,7 @@ 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: */ @@ -4960,5 +4960,5 @@ /*}}}*/ /*FUNCTION Tria::InputUpdateFromSolutionAdjointBalancethickness {{{*/ -void Tria::InputUpdateFromSolutionAdjointBalancethickness(double* solution){ +void Tria::InputUpdateFromSolutionAdjointBalancethickness(IssmDouble* solution){ const int numdof=NDOF1*NUMVERTICES; @@ -4966,6 +4966,6 @@ int i; int* doflist=NULL; - double values[numdof]; - double lambda[NUMVERTICES]; + IssmDouble values[numdof]; + IssmDouble lambda[NUMVERTICES]; /*Get dof list: */ @@ -5017,7 +5017,7 @@ /*}}}*/ /*FUNCTION Tria::SetControlInputsFromVector{{{*/ -void Tria::SetControlInputsFromVector(double* vector,int control_enum,int control_index){ - - double values[NUMVERTICES]; +void Tria::SetControlInputsFromVector(IssmDouble* vector,int control_enum,int control_index){ + + IssmDouble values[NUMVERTICES]; int doflist1[NUMVERTICES]; Input *input = NULL; @@ -5057,12 +5057,12 @@ /*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; @@ -5115,18 +5115,18 @@ /*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; @@ -5221,9 +5221,9 @@ /*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; @@ -5269,6 +5269,6 @@ int i; int* doflist=NULL; - double watercolumn; - double values[numdof]; + IssmDouble watercolumn; + IssmDouble values[numdof]; GaussTria* gauss=NULL; @@ -5299,5 +5299,5 @@ /*}}}*/ /*FUNCTION Tria::InputUpdateFromSolutionHydrology{{{*/ -void Tria::InputUpdateFromSolutionHydrology(double* solution){ +void Tria::InputUpdateFromSolutionHydrology(IssmDouble* solution){ /*Intermediaries*/ @@ -5306,5 +5306,5 @@ int i; int* doflist=NULL; - double values[numdof]; + IssmDouble values[numdof]; /*Get dof list: */ @@ -5315,5 +5315,5 @@ values[i]=solution[doflist[i]]; if(isnan(values[i])) _error_("NaN found in solution vector"); - if (values[i]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(); @@ -5441,13 +5441,13 @@ } /*}}}*/ -/*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*/ @@ -5466,9 +5466,9 @@ 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); @@ -5511,15 +5511,15 @@ 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; @@ -5620,10 +5620,10 @@ /*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; @@ -5688,7 +5688,7 @@ /*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; @@ -5731,7 +5731,7 @@ /*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; Index: /issm/trunk-jpl/src/c/objects/Elements/Tria.h =================================================================== --- /issm/trunk-jpl/src/c/objects/Elements/Tria.h (revision 12470) +++ /issm/trunk-jpl/src/c/objects/Elements/Tria.h (revision 12471) @@ -55,15 +55,15 @@ /*}}}*/ /*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); - #endif - void InputUpdateFromConstant(double constant, int name); + void InputUpdateFromMatrixDakota(IssmDouble* matrix, int nows, int ncols, int name, int type); + #endif + void InputUpdateFromConstant(IssmDouble constant, int name); void InputUpdateFromConstant(int constant, int name); void InputUpdateFromConstant(bool constant, int name); @@ -71,5 +71,5 @@ /*}}}*/ /*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); @@ -85,56 +85,56 @@ 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); @@ -148,17 +148,17 @@ 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 @@ -180,5 +180,5 @@ 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); @@ -186,14 +186,14 @@ 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_ @@ -207,6 +207,6 @@ 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 @@ -217,6 +217,6 @@ ElementVector* CreatePVectorAdjointStokes(void); ElementVector* CreatePVectorAdjointBalancethickness(void); - void InputUpdateFromSolutionAdjointBalancethickness( double* solution); - void InputUpdateFromSolutionAdjointHoriz( double* solution); + void InputUpdateFromSolutionAdjointBalancethickness( IssmDouble* solution); + void InputUpdateFromSolutionAdjointHoriz( IssmDouble* solution); #endif @@ -226,5 +226,5 @@ void CreateHydrologyWaterVelocityInput(void); void GetSolutionFromInputsHydrology(Vector* solution); - void InputUpdateFromSolutionHydrology(double* solution); + void InputUpdateFromSolutionHydrology(IssmDouble* solution); #endif #ifdef _HAVE_BALANCED_ Index: /issm/trunk-jpl/src/c/objects/Elements/TriaRef.cpp =================================================================== --- /issm/trunk-jpl/src/c/objects/Elements/TriaRef.cpp (revision 12470) +++ /issm/trunk-jpl/src/c/objects/Elements/TriaRef.cpp (revision 12471) @@ -57,5 +57,5 @@ /*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 @@ -70,5 +70,5 @@ int i; - double dbasis[NDOF2][NUMNODES]; + IssmDouble dbasis[NDOF2][NUMNODES]; /*Get dh1dh2dh3 in actual coordinate system: */ @@ -87,5 +87,5 @@ /*}}}*/ /*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. @@ -101,5 +101,5 @@ /*Same thing in the actual coordinate system: */ - double dbasis[NDOF2][NUMNODES]; + IssmDouble dbasis[NDOF2][NUMNODES]; /*Get dh1dh2dh3 in actual coordinates system : */ @@ -118,5 +118,5 @@ /*}}}*/ /*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] * @@ -126,5 +126,5 @@ */ - double l1l3[NUMNODES]; + IssmDouble l1l3[NUMNODES]; GetNodalFunctions(&l1l3[0],gauss); @@ -137,5 +137,5 @@ /*}}}*/ /*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] * @@ -145,5 +145,5 @@ */ - double l1l3[NUMNODES]; + IssmDouble l1l3[NUMNODES]; GetNodalFunctions(&l1l3[0],gauss); @@ -156,5 +156,5 @@ /*}}}*/ /*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 @@ -167,5 +167,5 @@ */ - double basis[NUMNODES]; + IssmDouble basis[NUMNODES]; /*Get dh1dh2dh3 in actual coordinate system: */ @@ -180,5 +180,5 @@ /*}}}*/ /*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. @@ -194,5 +194,5 @@ /*Same thing in the actual coordinate system: */ - double dbasis[NDOF2][NUMNODES]; + IssmDouble dbasis[NDOF2][NUMNODES]; /*Get dh1dh2dh3 in actual coordinates system : */ @@ -211,5 +211,5 @@ /*}}}*/ /*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. @@ -226,5 +226,5 @@ /*Same thing in the actual coordinate system: */ - double dbasis[NDOF2][NUMNODES]; + IssmDouble dbasis[NDOF2][NUMNODES]; /*Get dh1dh2dh3 in actual coordinates system : */ @@ -245,5 +245,5 @@ /*}}}*/ /*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 @@ -257,5 +257,5 @@ /*Same thing in the actual coordinate system: */ - double dbasis[NDOF2][NUMNODES]; + IssmDouble dbasis[NDOF2][NUMNODES]; /*Get dh1dh2dh3 in actual coordinates system : */ @@ -270,5 +270,5 @@ /*}}}*/ /*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 @@ -285,5 +285,5 @@ int i; - double basis[3]; + IssmDouble basis[3]; /*Get basis in actual coordinate system: */ @@ -307,8 +307,8 @@ /*}}}*/ /*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); @@ -327,8 +327,8 @@ /*}}}*/ /*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); @@ -343,8 +343,8 @@ /*}}}*/ /*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*/ @@ -358,9 +358,9 @@ /*}}}*/ /*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); @@ -380,8 +380,8 @@ /*}}}*/ /*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:*/ @@ -394,5 +394,5 @@ /*}}}*/ /*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.*/ @@ -404,8 +404,8 @@ /*}}}*/ /*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); @@ -418,11 +418,11 @@ /*}}}*/ /*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: */ @@ -445,5 +445,5 @@ /*}}}*/ /*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. */ @@ -464,5 +464,5 @@ /*}}}*/ /*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 @@ -475,5 +475,5 @@ /*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: */ @@ -487,5 +487,5 @@ /*}}}*/ /*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 @@ -493,5 +493,5 @@ /*nodal functions annd output: */ - double basis[3]; + IssmDouble basis[3]; /*Get nodal functions*/ Index: /issm/trunk-jpl/src/c/objects/Elements/TriaRef.h =================================================================== --- /issm/trunk-jpl/src/c/objects/Elements/TriaRef.h (revision 12470) +++ /issm/trunk-jpl/src/c/objects/Elements/TriaRef.h (revision 12471) @@ -25,24 +25,24 @@ /*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); };