Index: /issm/trunk/src/c/EnumDefinitions/EnumDefinitions.h =================================================================== --- /issm/trunk/src/c/EnumDefinitions/EnumDefinitions.h (revision 9270) +++ /issm/trunk/src/c/EnumDefinitions/EnumDefinitions.h (revision 9271) @@ -484,5 +484,5 @@ InputfilenameEnum, SpcvzEnum, - list, + Fake33Enum, NumberOfNodes2DEnum, NodeOnStokesEnum, @@ -548,5 +548,7 @@ ViscousHeatingEnum, QmuTemperatureEnum, - QmuRheologyBEnum + QmuRheologyBEnum, + HydrologyWaterVxEnum, + HydrologyWaterVyEnum }; Index: /issm/trunk/src/c/modules/EnumToStringx/EnumToStringx.cpp =================================================================== --- /issm/trunk/src/c/modules/EnumToStringx/EnumToStringx.cpp (revision 9270) +++ /issm/trunk/src/c/modules/EnumToStringx/EnumToStringx.cpp (revision 9271) @@ -426,4 +426,5 @@ case InputfilenameEnum : return "Inputfilename"; case SpcvzEnum : return "Spcvz"; + case Fake33Enum : return "Fake33"; case NumberOfNodes2DEnum : return "NumberOfNodes2D"; case NodeOnStokesEnum : return "NodeOnStokes"; @@ -490,4 +491,6 @@ case QmuTemperatureEnum : return "QmuTemperature"; case QmuRheologyBEnum : return "QmuRheologyB"; + case HydrologyWaterVxEnum : return "HydrologyWaterVx"; + case HydrologyWaterVyEnum : return "HydrologyWaterVy"; default : return "unknown"; Index: /issm/trunk/src/c/modules/ModelProcessorx/Hydrology/CreateConstraintsHydrology.cpp =================================================================== --- /issm/trunk/src/c/modules/ModelProcessorx/Hydrology/CreateConstraintsHydrology.cpp (revision 9270) +++ /issm/trunk/src/c/modules/ModelProcessorx/Hydrology/CreateConstraintsHydrology.cpp (revision 9271) @@ -4,4 +4,5 @@ #include "../../../Container/Container.h" +#include "../../../modules/modules.h" #include "../../../toolkits/toolkits.h" #include "../../../EnumDefinitions/EnumDefinitions.h" @@ -13,11 +14,6 @@ void CreateConstraintsHydrology(Constraints** pconstraints, IoModel* iomodel,FILE* iomodel_handle){ - /*Intermediary*/ - int i; - int count; - /*Output*/ Constraints *constraints = NULL; - Spc *spc = NULL; /*Recover pointer: */ @@ -26,27 +22,5 @@ /*Create constraints if they do not exist yet*/ if(!constraints) constraints = new Constraints(ConstraintsEnum); - - /*Fetch data: */ - IoModelFetchData(&iomodel->spcwatercolumn,NULL,NULL,iomodel_handle,SpcwatercolumnEnum); - - /*Initialize counter*/ - count=0; - - /*Create spcs from x,y,z, as well as the spc values on those spcs: */ - for (i=0;inumberofvertices;i++){ - - /*keep only this partition's nodes:*/ - if((iomodel->my_vertices[i])){ - - if ((int)iomodel->spcwatercolumn[2*i]){ - - constraints->AddObject(new Spc(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,1,*(iomodel->spcwatercolumn+2*i+1),HydrologyAnalysisEnum)); - count++; - } - } //if((my_vertices[i])) - } - - /*Free data: */ - xfree((void**)&iomodel->spcwatercolumn); + IoModelToConstraintsx(constraints,iomodel,iomodel_handle,SpcwatercolumnEnum,HydrologyAnalysisEnum); /*Assign output pointer: */ Index: /issm/trunk/src/c/modules/StringToEnumx/StringToEnumx.cpp =================================================================== --- /issm/trunk/src/c/modules/StringToEnumx/StringToEnumx.cpp (revision 9270) +++ /issm/trunk/src/c/modules/StringToEnumx/StringToEnumx.cpp (revision 9271) @@ -424,4 +424,5 @@ else if (strcmp(name,"Inputfilename")==0) return InputfilenameEnum; else if (strcmp(name,"Spcvz")==0) return SpcvzEnum; + else if (strcmp(name,"Fake33")==0) return Fake33Enum; else if (strcmp(name,"NumberOfNodes2D")==0) return NumberOfNodes2DEnum; else if (strcmp(name,"NodeOnStokes")==0) return NodeOnStokesEnum; @@ -488,4 +489,6 @@ else if (strcmp(name,"QmuTemperature")==0) return QmuTemperatureEnum; else if (strcmp(name,"QmuRheologyB")==0) return QmuRheologyBEnum; + else if (strcmp(name,"HydrologyWaterVx")==0) return HydrologyWaterVxEnum; + else if (strcmp(name,"HydrologyWaterVy")==0) return HydrologyWaterVyEnum; else _error_("Enum %s not found",name); Index: /issm/trunk/src/c/objects/Elements/Tria.cpp =================================================================== --- /issm/trunk/src/c/objects/Elements/Tria.cpp (revision 9270) +++ /issm/trunk/src/c/objects/Elements/Tria.cpp (revision 9271) @@ -376,4 +376,50 @@ } /*}}}*/ +/*FUNCTION Tria::CreateHydrologyWaterVelocityInput {{{1*/ +void Tria::CreateHydrologyWaterVelocityInput(void){ + + /*material parameters: */ + double mu_water=0.001787; //unit= [N s/m2] + double w; + double rho_ice, rho_water, g; + double dsdx,dsdy,dbdx,dbdy; + double vx[NUMVERTICES]; + double vy[NUMVERTICES]; + GaussTria *gauss = NULL; + + /*Retrieve all inputs and parameters*/ + rho_ice=matpar->GetRhoIce(); + rho_water=matpar->GetRhoWater(); + g=matpar->GetG(); + Input* surfaceslopex_input=inputs->GetInput(SurfaceSlopeXEnum); _assert_(surfaceslopex_input); + Input* surfaceslopey_input=inputs->GetInput(SurfaceSlopeYEnum); _assert_(surfaceslopey_input); + Input* bedslopex_input=inputs->GetInput(BedSlopeXEnum); _assert_(bedslopex_input); + Input* bedslopey_input=inputs->GetInput(BedSlopeYEnum); _assert_(bedslopey_input); + Input* watercolumn_input=inputs->GetInput(WatercolumnEnum); _assert_(watercolumn_input); + + gauss=new GaussTria(); + for (int iv=0;ivGaussVertex(iv); + surfaceslopex_input->GetParameterValue(&dsdx,gauss); + surfaceslopey_input->GetParameterValue(&dsdy,gauss); + bedslopex_input->GetParameterValue(&dbdx,gauss); + bedslopey_input->GetParameterValue(&dbdy,gauss); + watercolumn_input->GetParameterValue(&w,gauss); + + /* Water velocity x and y components */ + vx[iv]= pow(w,2)/(12 * mu_water)*(rho_ice*g*dsdx+(rho_water-rho_ice)*g*dbdx); + vy[iv]= pow(w,2)/(12 * mu_water)*(rho_ice*g*dsdy+(rho_water-rho_ice)*g*dbdy); + //vx[iv]=0.001; + //vy[iv]=0; + } + + /*clean-up*/ + delete gauss; + + /*Add to inputs*/ + this->inputs->AddInput(new TriaVertexInput(HydrologyWaterVxEnum,vx)); + this->inputs->AddInput(new TriaVertexInput(HydrologyWaterVyEnum,vy)); +} +/*}}}*/ /*FUNCTION Tria::CreateKMatrix {{{1*/ void Tria::CreateKMatrix(Mat Kff, Mat Kfs,Vec df){ @@ -783,4 +829,107 @@ } /*}}}*/ +/*FUNCTION Tria::CreateKMatrixHydrology{{{1*/ +ElementMatrix* Tria::CreateKMatrixHydrology(void){ + + /*Constants*/ + const int numdof=NDOF1*NUMVERTICES; + + /*Intermediaries */ + int artdiff; + 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}; + GaussTria *gauss=NULL; + + /*Initialize Element matrix*/ + ElementMatrix* Ke=new ElementMatrix(nodes,NUMVERTICES,this->parameters,NoneApproximationEnum); + + /*Create water velocity vx and vy from current inputs*/ + CreateHydrologyWaterVelocityInput(); + + /*Retrieve all inputs and parameters*/ + GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES); + this->parameters->FindParam(&dt,DtEnum); + this->parameters->FindParam(&artdiff,ArtDiffEnum); + Input* vx_input=inputs->GetInput(HydrologyWaterVxEnum); _assert_(vx_input); + Input* vy_input=inputs->GetInput(HydrologyWaterVyEnum); _assert_(vy_input); + h=this->GetArea(); + + /* Start looping on the number of gaussian points: */ + gauss=new GaussTria(2); + for (ig=gauss->begin();igend();ig++){ + + gauss->GaussPoint(ig); + + GetJacobianDeterminant2d(&Jdettria, &xyz_list[0][0],gauss); + GetL(&L[0], &xyz_list[0][0], gauss,NDOF1); + + vx_input->GetParameterValue(&vx,gauss); + vy_input->GetParameterValue(&vy,gauss); + vx_input->GetParameterDerivativeValue(&dvx[0],&xyz_list[0][0],gauss); + vy_input->GetParameterDerivativeValue(&dvy[0],&xyz_list[0][0],gauss); + + DL_scalar=gauss->weight*Jdettria; + + TripleMultiply( &L[0],1,numdof,1, + &DL_scalar,1,1,0, + &L[0],1,numdof,0, + &Ke->values[0],1); + + GetBPrognostic(&B[0][0], &xyz_list[0][0], gauss); + GetBprimePrognostic(&Bprime[0][0], &xyz_list[0][0], gauss); + + dvxdx=dvx[0]; + dvydy=dvy[1]; + DL_scalar=dt*gauss->weight*Jdettria; + + DL[0][0]=DL_scalar*dvxdx; + DL[1][1]=DL_scalar*dvydy; + DLprime[0][0]=DL_scalar*vx; + DLprime[1][1]=DL_scalar*vy; + + TripleMultiply( &B[0][0],2,numdof,1, + &DL[0][0],2,2,0, + &B[0][0],2,numdof,0, + &Ke->values[0],1); + + TripleMultiply( &B[0][0],2,numdof,1, + &DLprime[0][0],2,2,0, + &Bprime[0][0],2,numdof,0, + &Ke->values[0],1); + + if(artdiff){ + vel=sqrt(pow(vx,2.)+pow(vy,2.)); + K[0][0]=h/(2*vel)*vx*vx; + K[1][0]=h/(2*vel)*vy*vx; + K[0][1]=h/(2*vel)*vx*vy; + K[1][1]=h/(2*vel)*vy*vy; + KDL[0][0]=DL_scalar*K[0][0]; + KDL[1][0]=DL_scalar*K[1][0]; + KDL[0][1]=DL_scalar*K[0][1]; + KDL[1][1]=DL_scalar*K[1][1]; + + TripleMultiply( &Bprime[0][0],2,numdof,1, + &KDL[0][0],2,2,0, + &Bprime[0][0],2,numdof,0, + &Ke->values[0],1); + } + } + + /*Clean up and return*/ + delete gauss; + return Ke; +} +/*}}}*/ /*FUNCTION Tria::CreateKMatrixMelting {{{1*/ ElementMatrix* Tria::CreateKMatrixMelting(void){ @@ -820,83 +969,4 @@ &L[0],1,numdof,0, &Ke->values[0],1); - } - - /*Clean up and return*/ - delete gauss; - return Ke; -} -/*}}}*/ -/*FUNCTION Tria::CreateKMatrixHydrology {{{1*/ -ElementMatrix* Tria::CreateKMatrixHydrology(void){ - - /*Constants*/ - const int numdof=NDOF1*NUMVERTICES; - - /*Intermediaries */ - int i,j,ig; - double Jdettria,dt; - double xyz_list[NUMVERTICES][3]; - double L[NUMVERTICES]; - double dL[2][NUMVERTICES]; - double dLx[NUMVERTICES]; - double dLy[NUMVERTICES]; - double K[2]; - double hydro_p,hydro_q,hydro_CR,hydro_n,rho_water,g,gamma; - double DL_scalar1; - double DL_scalar2; - double DL_scalar3; - GaussTria *gauss=NULL; - - /*Initialize Element matrix*/ - ElementMatrix* Ke=new ElementMatrix(nodes,NUMVERTICES,this->parameters,NoneApproximationEnum); - - /*Retrieve all inputs and parameters*/ - GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES); - this->parameters->FindParam(&dt,DtEnum); - - /*retrieve material parameters: */ - g =matpar->GetG(); - rho_water=matpar->GetRhoWater(); - hydro_p =matpar->GetHydroP(); - hydro_q =matpar->GetHydroQ(); - hydro_CR =matpar->GetHydroCR(); - hydro_n =matpar->GetHydroN(); - gamma=1/(hydro_n*pow(hydro_CR,hydro_p)*pow((rho_water*g),hydro_q)); - - /* Start looping on the number of gaussian points: */ - gauss=new GaussTria(2); - for (ig=gauss->begin();igend();ig++){ - - gauss->GaussPoint(ig); - - GetJacobianDeterminant2d(&Jdettria, &xyz_list[0][0],gauss); - GetNodalFunctions(&L[0],gauss); - GetNodalFunctionsDerivatives(&dL[0][0],&xyz_list[0][0],gauss); - for(i=0;iweight*Jdettria; - TripleMultiply( &L[0],1,numdof,1, - &DL_scalar1,1,1,0, - &L[0],1,numdof,0, - &Ke->values[0],1); - } - - if(dt){ - DL_scalar2=-gauss->weight*Jdettria*gamma*K[0]*dt; - DL_scalar3=-gauss->weight*Jdettria*gamma*K[1]*dt; - - } - else{ - DL_scalar2=-gauss->weight*Jdettria*gamma*K[0]; - DL_scalar3=-gauss->weight*Jdettria*gamma*K[1]; - } - - TripleMultiply(&dLx[0],1,numdof,1, &DL_scalar2,1,1,0, &L[0],1,numdof,0,&Ke->values[0],1); - TripleMultiply(&dLy[0],1,numdof,1, &DL_scalar3,1,1,0, &L[0],1,numdof,0,&Ke->values[0],1); } @@ -1865,5 +1935,5 @@ double old_watercolumn_g; double xyz_list[NUMVERTICES][3]; - double L[NUMVERTICES]; + double basis[numdof]; GaussTria* gauss=NULL; @@ -1874,5 +1944,5 @@ GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES); this->parameters->FindParam(&dt,DtEnum); - Input* basal_melting_input=inputs->GetInput(BasalMeltingRateEnum); _assert_(basal_melting_input); + Input* basal_melting_input=inputs->GetInput(BasalMeltingRateEnum);_assert_(basal_melting_input); Input* old_watercolumn_input=inputs->GetInput(WaterColumnOldEnum);_assert_(old_watercolumn_input); @@ -1885,13 +1955,11 @@ GetJacobianDeterminant2d(&Jdettria, &xyz_list[0][0],gauss); - GetL(&L[0], &xyz_list[0][0], gauss,NDOF1); + GetNodalFunctions(basis, gauss); basal_melting_input->GetParameterValue(&basal_melting_g,gauss); old_watercolumn_input->GetParameterValue(&old_watercolumn_g,gauss); - //printf("ig=%d, old_watercolumn_g=%f,Jdettria=%f\n",ig,old_watercolumn_g,Jdettria);//bk - //if(old_watercolumn_g<0) printf("old_watercolumn_g=%g,Jdettria=%f\n",old_watercolumn_g,Jdettria);//bk - - if(dt)for(i=0;ivalues[i]+=Jdettria*gauss->weight*(old_watercolumn_g+dt*basal_melting_g)*L[i]; - else for(i=0;ivalues[i]+=Jdettria*gauss->weight*basal_melting_g*L[i]; + + if(dt)for(i=0;ivalues[i]+=Jdettria*gauss->weight*(old_watercolumn_g+dt*basal_melting_g)*basis[i]; + else for(i=0;ivalues[i]+=Jdettria*gauss->weight*basal_melting_g*basis[i]; } @@ -2334,76 +2402,4 @@ return &this->horizontalneighborsids[0]; -} -/*}}}*/ -/*FUNCTION Tria::GetHydrologyK {{{1*/ -void Tria::GetHydrologyK(double* K,double* xyz_list,GaussTria* gauss){ - - /*material parameters: */ - double rho_ice, rho_water, g,hydro_p,hydro_q; - double dsdx, dsdy; - double dbdx, dbdy; - double surface_slope; - double kn, w; - - /*Retrieve all inputs and parameters*/ - rho_ice=matpar->GetRhoIce(); - rho_water=matpar->GetRhoWater(); - g=matpar->GetG(); - kn=matpar->GetKn(); - hydro_p=matpar->GetHydroP(); - hydro_q=matpar->GetHydroQ(); - Input* surfaceslopex_input=inputs->GetInput(SurfaceSlopeXEnum); _assert_(surfaceslopex_input); - Input* surfaceslopey_input=inputs->GetInput(SurfaceSlopeYEnum); _assert_(surfaceslopey_input); - Input* bedslopex_input=inputs->GetInput(BedSlopeXEnum); _assert_(bedslopex_input); - Input* bedslopey_input=inputs->GetInput(BedSlopeYEnum); _assert_(bedslopey_input); - Input* watercolumn_input=inputs->GetInput(WatercolumnEnum); _assert_(watercolumn_input); - - /*recover slopes: */ - surfaceslopex_input->GetParameterValue(&dsdx,gauss); - surfaceslopey_input->GetParameterValue(&dsdy,gauss); - bedslopex_input->GetParameterValue(&dbdx,gauss); - bedslopey_input->GetParameterValue(&dbdy,gauss); - - /*magnitude of surface slope: */ - surface_slope=sqrt( pow(dsdx,2) + pow(dsdy,2)); - - /*recover water column: */ - watercolumn_input->GetParameterValue(&w,gauss); - - K[0]=pow(w,hydro_p)*pow((rho_ice*g*dsdx+(rho_water/rho_ice-1)*rho_ice*g*dbdx) - rho_ice * g * kn/w * (dsdx - dbdx ) * surface_slope,hydro_q); - K[1]=pow(w,hydro_p)*pow((rho_ice*g*dsdy+(rho_water/rho_ice-1)*rho_ice*g*dbdy) - rho_ice * g * kn/w * (dsdy - dbdy ) * surface_slope,hydro_q); - -} -/*}}}*/ -/*FUNCTION Tria::GetHydrologyWaterVelocity {{{1*/ -void Tria::GetHydrologyWaterVelocity(double* v,GaussTria* gauss){ - /*material parameters: */ - double mu_water=0.001787; //unit= [N s/m2] - double w; - double rho_ice, rho_water, g; - double dsdx,dsdy,dbdx,dbdy; - - /*Retrieve all inputs and parameters*/ - rho_ice=matpar->GetRhoIce(); - rho_water=matpar->GetRhoWater(); - g=matpar->GetG(); - Input* surfaceslopex_input=inputs->GetInput(SurfaceSlopeXEnum); _assert_(surfaceslopex_input); - Input* surfaceslopey_input=inputs->GetInput(SurfaceSlopeYEnum); _assert_(surfaceslopey_input); - Input* bedslopex_input=inputs->GetInput(BedSlopeXEnum); _assert_(bedslopex_input); - Input* bedslopey_input=inputs->GetInput(BedSlopeYEnum); _assert_(bedslopey_input); - Input* watercolumn_input=inputs->GetInput(WatercolumnEnum); _assert_(watercolumn_input); - - /*recover slopes: */ - surfaceslopex_input->GetParameterValue(&dsdx,gauss); - surfaceslopey_input->GetParameterValue(&dsdy,gauss); - bedslopex_input->GetParameterValue(&dbdx,gauss); - bedslopey_input->GetParameterValue(&dbdy,gauss); - - /*recover water column: */ - watercolumn_input->GetParameterValue(&w,gauss); - - /* Water velocity x and y components */ - v[0]=pow(w,2)/(12 * mu_water)*(rho_ice*g*dsdx+(rho_water-rho_ice)*g*dbdx); - v[1]=pow(w,2)/(12 * mu_water)*(rho_ice*g*dsdy+(rho_water-rho_ice)*g*dbdy); } /*}}}*/ Index: /issm/trunk/src/c/objects/Elements/Tria.h =================================================================== --- /issm/trunk/src/c/objects/Elements/Tria.h (revision 9270) +++ /issm/trunk/src/c/objects/Elements/Tria.h (revision 9271) @@ -181,5 +181,5 @@ int GetElementType(void); void GetDofList(int** pdoflist,int approximation_enum,int setenum); - void GetHydrologyWaterVelocity(double* v, GaussTria* gauss); + void CreateHydrologyWaterVelocityInput(void); void GetDofList1(int* doflist); void GetSidList(int* sidlist); @@ -202,5 +202,4 @@ void SetClone(int* minranks); void SurfaceNormal(double* surface_normal, double xyz_list[3][3]); - void GetHydrologyK(double* K,double* xyz_list,GaussTria* gauss); /*}}}*/ Index: /issm/trunk/src/c/objects/IoModel.cpp =================================================================== --- /issm/trunk/src/c/objects/IoModel.cpp (revision 9270) +++ /issm/trunk/src/c/objects/IoModel.cpp (revision 9271) @@ -223,11 +223,4 @@ IoModelFetchData(&this->gl_melting_rate,iomodel_handle,GlMeltingRateEnum); IoModelFetchData(&this->rheology_law,iomodel_handle,RheologyLawEnum); - - /*!Get hydrology parameters: */ - IoModelFetchData(&this->hydro_kn,iomodel_handle,HydroKnEnum); - IoModelFetchData(&this->hydro_p,iomodel_handle,HydroPEnum); - IoModelFetchData(&this->hydro_q,iomodel_handle,HydroQEnum); - IoModelFetchData(&this->hydro_CR,iomodel_handle,HydroCREnum); - IoModelFetchData(&this->hydro_n,iomodel_handle,HydroNEnum); /*qmu: */ Index: /issm/trunk/src/c/objects/IoModel.h =================================================================== --- /issm/trunk/src/c/objects/IoModel.h (revision 9270) +++ /issm/trunk/src/c/objects/IoModel.h (revision 9271) @@ -215,11 +215,4 @@ int loadcounter; //keep track of how many loads are being created in each analysis type int constraintcounter; //keep track of how many constraints are being created in each analysis type - - /*hydrology parameters: */ - double hydro_kn; - double hydro_p; - double hydro_q; - double hydro_CR; - double hydro_n; /*}}}*/ /*Methods: {{{1*/ Index: /issm/trunk/src/c/objects/Materials/Matpar.cpp =================================================================== --- /issm/trunk/src/c/objects/Materials/Matpar.cpp (revision 9270) +++ /issm/trunk/src/c/objects/Materials/Matpar.cpp (revision 9271) @@ -37,9 +37,4 @@ this->thermal_exchange_velocity = iomodel->thermal_exchange_velocity; this->g = iomodel->g; - this->kn = iomodel->hydro_kn; - this->hydro_p = iomodel->hydro_p; - this->hydro_q = iomodel->hydro_q; - this->hydro_CR = iomodel->hydro_CR; - this->hydro_n = iomodel->hydro_n; } /*}}}1*/ @@ -361,29 +356,4 @@ } /*}}}1*/ -/*FUNCTION Matpar::GetKn {{{1*/ -double Matpar::GetKn(){ - return kn; -} -/*}}}1*/ -/*FUNCTION Matpar::GetHydroP {{{1*/ -double Matpar::GetHydroP(){ - return hydro_p; -} -/*}}}1*/ -/*FUNCTION Matqar::GetHydroQ {{{1*/ -double Matpar::GetHydroQ(){ - return hydro_q; -} -/*}}}1*/ -/*FUNCTION Matpar::GetHydroCR {{{1*/ -double Matpar::GetHydroCR(){ - return hydro_CR; -} -/*}}}1*/ -/*FUNCTION Matpar::GetHydroN {{{1*/ -double Matpar::GetHydroN(){ - return hydro_n; -} -/*}}}1*/ /*FUNCTION Matpar::TMeltingPoint {{{1*/ double Matpar::TMeltingPoint(double pressure){ Index: /issm/trunk/src/c/objects/Materials/Matpar.h =================================================================== --- /issm/trunk/src/c/objects/Materials/Matpar.h (revision 9270) +++ /issm/trunk/src/c/objects/Materials/Matpar.h (revision 9271) @@ -27,11 +27,4 @@ double thermal_exchange_velocity; double g; - - /*hydrology: */ - double kn; - double hydro_p; - double hydro_q; - double hydro_CR; - double hydro_n; public: Index: /issm/trunk/src/c/solutions/hydrology_core.cpp =================================================================== --- /issm/trunk/src/c/solutions/hydrology_core.cpp (revision 9270) +++ /issm/trunk/src/c/solutions/hydrology_core.cpp (revision 9271) @@ -45,4 +45,5 @@ if(nsteps)_printf_(VerboseSolution(),"time step:%i/%i\n",i+1,nsteps); time=(i+1)*dt; + femmodel->parameters->SetParam(time,TimeEnum); /*call hydrology_core step: */ @@ -52,4 +53,6 @@ _printf_(VerboseSolution()," saving results\n"); InputToResultx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,WatercolumnEnum,i+1,time); + InputToResultx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,HydrologyWaterVxEnum,i+1,time); + InputToResultx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,HydrologyWaterVyEnum,i+1,time); } Index: /issm/trunk/src/c/solvers/solver_adjoint_linear.cpp =================================================================== --- /issm/trunk/src/c/solvers/solver_adjoint_linear.cpp (revision 9270) +++ /issm/trunk/src/c/solvers/solver_adjoint_linear.cpp (revision 9271) @@ -22,4 +22,5 @@ /*Recover parameters: */ femmodel->parameters->FindParam(&configuration_type,ConfigurationTypeEnum); + UpdateConstraintsx(femmodel->nodes,femmodel->constraints,femmodel->parameters); SystemMatricesx(&Kff, &Kfs, &pf, &df, NULL,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters); Index: /issm/trunk/src/c/solvers/solver_nonlinear.cpp =================================================================== --- /issm/trunk/src/c/solvers/solver_nonlinear.cpp (revision 9270) +++ /issm/trunk/src/c/solvers/solver_nonlinear.cpp (revision 9271) @@ -31,9 +31,9 @@ int configuration_type; - /*Recover parameters: */ femmodel->parameters->FindParam(&min_mechanical_constraints,MinMechanicalConstraintsEnum); femmodel->parameters->FindParam(&max_nonlinear_iterations,MaxNonlinearIterationsEnum); femmodel->parameters->FindParam(&configuration_type,ConfigurationTypeEnum); + UpdateConstraintsx(femmodel->nodes,femmodel->constraints,femmodel->parameters); /*Were loads requested as output? : */ Index: /issm/trunk/src/c/solvers/solver_stokescoupling_nonlinear.cpp =================================================================== --- /issm/trunk/src/c/solvers/solver_stokescoupling_nonlinear.cpp (revision 9270) +++ /issm/trunk/src/c/solvers/solver_stokescoupling_nonlinear.cpp (revision 9271) @@ -35,4 +35,5 @@ femmodel->parameters->FindParam(&min_mechanical_constraints,MinMechanicalConstraintsEnum); femmodel->parameters->FindParam(&max_nonlinear_iterations,MaxNonlinearIterationsEnum); + UpdateConstraintsx(femmodel->nodes,femmodel->constraints,femmodel->parameters); count=1; Index: /issm/trunk/src/m/classes/model.m =================================================================== --- /issm/trunk/src/m/classes/model.m (revision 9270) +++ /issm/trunk/src/m/classes/model.m (revision 9271) @@ -268,11 +268,5 @@ basal_melting_rate_correction_apply = {0,true,'Integer'}; pressure = {NaN,true,'DoubleMat',1}; - %Hydrology watercolumn = {NaN,true,'DoubleMat',1}; - hydro_n = {0,true,'Double'}; - hydro_CR = {0,true,'Double'}; - hydro_p = {0,true,'Double'}; - hydro_q = {0,true,'Double'}; - hydro_kn = {0,true,'Double'}; %Parallelisation @@ -849,12 +843,4 @@ md.petscoptions=addoptions(md.petscoptions,DiagnosticVertAnalysisEnum,mumpsoptions); - - %hydrology: from Johnson's 2002 thesis, section 3.5.4 - md.hydro_n=.02; - md.hydro_CR=2; - md.hydro_p=2; - md.hydro_q=1; - md.hydro_kn=0; - %Rheology law: what is the temperature dependence of B with T %available: None, Paterson and Arrhenius Index: /issm/trunk/src/m/enum/Fake33Enum.m =================================================================== --- /issm/trunk/src/m/enum/Fake33Enum.m (revision 9271) +++ /issm/trunk/src/m/enum/Fake33Enum.m (revision 9271) @@ -0,0 +1,11 @@ +function macro=Fake33Enum() +%FAKE33ENUM - Enum of Fake33 +% +% WARNING: DO NOT MODIFY THIS FILE +% this file has been automatically generated by src/c/EnumDefinitions/Synchronize.sh +% Please read src/c/EnumDefinitions/README for more information +% +% Usage: +% macro=Fake33Enum() + +macro=StringToEnum('Fake33'); Index: /issm/trunk/src/m/enum/HydrologyWaterVxEnum.m =================================================================== --- /issm/trunk/src/m/enum/HydrologyWaterVxEnum.m (revision 9271) +++ /issm/trunk/src/m/enum/HydrologyWaterVxEnum.m (revision 9271) @@ -0,0 +1,11 @@ +function macro=HydrologyWaterVxEnum() +%HYDROLOGYWATERVXENUM - Enum of HydrologyWaterVx +% +% WARNING: DO NOT MODIFY THIS FILE +% this file has been automatically generated by src/c/EnumDefinitions/Synchronize.sh +% Please read src/c/EnumDefinitions/README for more information +% +% Usage: +% macro=HydrologyWaterVxEnum() + +macro=StringToEnum('HydrologyWaterVx'); Index: /issm/trunk/src/m/enum/HydrologyWaterVyEnum.m =================================================================== --- /issm/trunk/src/m/enum/HydrologyWaterVyEnum.m (revision 9271) +++ /issm/trunk/src/m/enum/HydrologyWaterVyEnum.m (revision 9271) @@ -0,0 +1,11 @@ +function macro=HydrologyWaterVyEnum() +%HYDROLOGYWATERVYENUM - Enum of HydrologyWaterVy +% +% WARNING: DO NOT MODIFY THIS FILE +% this file has been automatically generated by src/c/EnumDefinitions/Synchronize.sh +% Please read src/c/EnumDefinitions/README for more information +% +% Usage: +% macro=HydrologyWaterVyEnum() + +macro=StringToEnum('HydrologyWaterVy');