Index: /issm/trunk-jpl/src/c/classes/FemModel.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/FemModel.cpp (revision 14966) +++ /issm/trunk-jpl/src/c/classes/FemModel.cpp (revision 14967) @@ -1459,5 +1459,5 @@ #ifdef _HAVE_DAKOTA_ void FemModel::DakotaResponsesx(double* d_responses,char** responses_descriptors,int numresponsedescriptors,int d_numresponses){/*{{{*/ - + int i,j; int my_rank; @@ -1562,9 +1562,9 @@ void FemModel::Deflection(Vector* wg,Vector* dwgdt, IssmDouble* x, IssmDouble* y){ /*{{{*/ - int i; - + int i; + /*intermediary: */ Element *element = NULL; - + /*Go through elements, and add contribution from each element to the deflection vector wg:*/ for (i=0;iSize();i++){ @@ -1575,2 +1575,22 @@ /*}}}*/ #endif + +void FemModel::BasisIntegralsx(void){ /*{{{*/ + + Vector* basisg=NULL; + + /*Vector allocation*/ + basisg=new Vector(nodes->NumberOfNodes()); + + for (int i=0;iSize();i++){ + Element* element=dynamic_cast(elements->GetObjectByOffset(i)); + element->BasisIntegral(basisg); + } + /*Assemble*/ + basisg->Assemble(); + + /*Update Inputs*/ + InputUpdateFromVectorx(elements,nodes,vertices,loads,materials,parameters,basisg,BasisIntegralEnum,NodesEnum); + +} +/*}}}*/ Index: /issm/trunk-jpl/src/c/classes/FemModel.h =================================================================== --- /issm/trunk-jpl/src/c/classes/FemModel.h (revision 14966) +++ /issm/trunk-jpl/src/c/classes/FemModel.h (revision 14967) @@ -100,4 +100,5 @@ void UpdateConstraintsx(void); int UpdateVertexPositionsx(void); + void BasisIntegralsx(void); }; Index: /issm/trunk-jpl/src/c/classes/objects/Elements/Element.h =================================================================== --- /issm/trunk-jpl/src/c/classes/objects/Elements/Element.h (revision 14966) +++ /issm/trunk-jpl/src/c/classes/objects/Elements/Element.h (revision 14967) @@ -128,5 +128,6 @@ #ifdef _HAVE_HYDROLOGY_ - virtual void GetHydrologyDCInefficientHmax(IssmDouble* ph_max, Node* innode)=0; + virtual void GetHydrologyDCInefficientHmax(IssmDouble* ph_max, Node* innode)=0; + virtual void BasisIntegral(Vector* basisg)=0; #endif }; Index: /issm/trunk-jpl/src/c/classes/objects/Elements/Penta.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/objects/Elements/Penta.cpp (revision 14966) +++ /issm/trunk-jpl/src/c/classes/objects/Elements/Penta.cpp (revision 14967) @@ -9347,3 +9347,9 @@ } /*}}}*/ +/*FUNCTION Tria::BasisIntegral {{{*/ +void Penta::BasisIntegral(Vector* basisg){ + _error_("Hydrological stuff not suported in Penta"); +} +/*}}}*/ + #endif Index: /issm/trunk-jpl/src/c/classes/objects/Elements/Penta.h =================================================================== --- /issm/trunk-jpl/src/c/classes/objects/Elements/Penta.h (revision 14966) +++ /issm/trunk-jpl/src/c/classes/objects/Elements/Penta.h (revision 14967) @@ -305,4 +305,5 @@ void CreateHydrologyWaterVelocityInput(void); void GetHydrologyDCInefficientHmax(IssmDouble* ph_max, Node* innode); + void BasisIntegral(Vector* gbasis); #endif #ifdef _HAVE_THERMAL_ Index: /issm/trunk-jpl/src/c/classes/objects/Elements/Tria.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/objects/Elements/Tria.cpp (revision 14966) +++ /issm/trunk-jpl/src/c/classes/objects/Elements/Tria.cpp (revision 14967) @@ -1881,29 +1881,43 @@ void Tria::InputUpdateFromVector(IssmDouble* vector, int name, int type){ + + const int numdof = NDOF1 *NUMVERTICES; + int *doflist = NULL; + IssmDouble values[numdof]; + /*Check that name is an element input*/ if (!IsInput(name)) return; switch(type){ - - case VertexEnum: { - - /*New TriaP1Input*/ - IssmDouble values[3]; - - /*Get values on the 3 vertices*/ - for (int i=0;i<3;i++){ - values[i]=vector[this->nodes[i]->GetVertexPid()]; - } - - /*update input*/ - if (name==MaterialsRheologyBbarEnum || name==MaterialsRheologyBEnum || name==MaterialsRheologyZEnum || name==MaterialsRheologyZbarEnum){ - material->inputs->AddInput(new TriaP1Input(name,values)); - } - else{ - this->inputs->AddInput(new TriaP1Input(name,values)); - } - return; + case VertexEnum: + /*Get values on the 3 vertices*/ + for (int i=0;i<3;i++){ + values[i]=vector[this->nodes[i]->GetVertexPid()]; } - default: + /*update input*/ + if (name==MaterialsRheologyBbarEnum || name==MaterialsRheologyBEnum || name==MaterialsRheologyZEnum || name==MaterialsRheologyZbarEnum){ + material->inputs->AddInput(new TriaP1Input(name,values)); + } + else{ + this->inputs->AddInput(new TriaP1Input(name,values)); + } + return; + + case NodesEnum: + /*Get dof list: */ + GetDofList(&doflist,NoneApproximationEnum,GsetEnum); + + /*Use the dof list to index into the vector: */ + for(int i=0;i(values[i])) _error_("NaN found in solution vector"); + } + /*Add input to the element: */ + this->inputs->AddInput(new TriaP1Input(BasisIntegralEnum,values)); + + /*Free ressources:*/ + xDelete(doflist); + default: + _error_("type " << type << " (" << EnumToStringx(type) << ") not implemented yet"); } @@ -2053,4 +2067,5 @@ name==OldGradientEnum || name==ConvergedEnum || + name==BasisIntegralEnum || name==QmuVxEnum || name==QmuVyEnum || @@ -6175,4 +6190,5 @@ /*Loading term*/ water_input->GetInputValue(&water_load,gauss); + if(this->id==1) printf("Qin %e \n ", water_load); scalar = Jdet*gauss->weight*water_load; if(reCast(dt)) scalar = scalar*dt; @@ -6232,4 +6248,5 @@ /*Loading term*/ residual_input->GetInputValue(&residual_load,gauss); + if(this->id==1) printf("Qres %e \n ", residual_load); scalar = Jdet*gauss->weight*residual_load; if(reCast(dt)) scalar = scalar*dt; @@ -6316,9 +6333,10 @@ const int numdof = NDOF1 *NUMVERTICES; int *doflist = NULL; - bool converged; + bool converged; IssmDouble values[numdof]; IssmDouble residual[numdof]; + IssmDouble intbasis[numdof]; IssmDouble sediment_storing; - IssmDouble penalty_factor; + IssmDouble penalty_factor, dt; IssmDouble kmax, kappa, h_max; @@ -6335,5 +6353,10 @@ /*If converged keep the residual in mind*/ this->inputs->GetInputValue(&converged,ConvergedEnum); + GetInputListOnVertices(&intbasis[0],BasisIntegralEnum); + + if(this->id==1) printf("val %e \n ", values[1]); + if(converged){ + this->parameters->FindParam(&dt,TimesteppingTimeStepEnum); this->parameters->FindParam(&kmax,HydrologySedimentKmaxEnum); this->parameters->FindParam(&penalty_factor,HydrologydcPenaltyFactorEnum); @@ -6344,9 +6367,13 @@ this->GetHydrologyDCInefficientHmax(&h_max,nodes[i]); if(values[i]>h_max) - residual[i]=kappa*(values[i]-h_max)*sediment_storing; + residual[i]=kappa*(values[i]-h_max)*sediment_storing/(dt*intbasis[i]); else residual[i]=0.0; } - } + if(this->id==1){ + printf("res %e val %e h-max %e Stor %e \n ", residual[1], values[1], h_max, sediment_storing); + printf("kappa %e kmax %e pen %e dt %e \n", kappa, kmax,penalty_factor,dt); + } + } /*Add input to the element: */ @@ -6421,4 +6448,40 @@ } /*}}}*/ +/*FUNCTION Tria::BasisIntegral {{{*/ +void Tria::BasisIntegral(Vector* basisg){ + + /*Constants*/ + const int numdof=NDOF1*NUMVERTICES; + + /*Intermediaries */ + IssmDouble Jdet; + IssmDouble xyz_list[NUMVERTICES][3]; + IssmDouble basis[numdof]; + IssmDouble basisint[numdof]={0.}; + int *doflist=NULL; + GaussTria* gauss=NULL; + + GetDofList(&doflist,NoneApproximationEnum,GsetEnum); + + /* Start looping on the number of gaussian points: */ + gauss=new GaussTria(2); + for(int ig=gauss->begin();igend();ig++){ + + gauss->GaussPoint(ig); + + GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss); + GetNodalFunctions(&basis[0], gauss); + + for(int i=0;iweight*basis[i]; + } + + basisg->SetValues(numdof,doflist,&basisint[0],ADD_VAL); + + /*Clean up and return*/ + delete gauss; + xDelete(doflist); +} +/*}}}*/ + #endif Index: /issm/trunk-jpl/src/c/classes/objects/Elements/Tria.h =================================================================== --- /issm/trunk-jpl/src/c/classes/objects/Elements/Tria.h (revision 14966) +++ /issm/trunk-jpl/src/c/classes/objects/Elements/Tria.h (revision 14967) @@ -253,4 +253,5 @@ void InputUpdateFromSolutionHydrologyDCEfficient(IssmDouble* solution); void GetHydrologyDCInefficientHmax(IssmDouble* ph_max, Node* innode); + void BasisIntegral(Vector* gbasis); #endif #ifdef _HAVE_BALANCED_ Index: /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h (revision 14966) +++ /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h (revision 14967) @@ -108,4 +108,5 @@ HydrologyEfficientEnum, HydrologySedimentKmaxEnum, + BasisIntegralEnum, IndependentObjectEnum, InversionControlParametersEnum, Index: /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp (revision 14966) +++ /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp (revision 14967) @@ -116,4 +116,5 @@ case HydrologyEfficientEnum : return "HydrologyEfficient"; case HydrologySedimentKmaxEnum : return "HydrologySedimentKmax"; + case BasisIntegralEnum : return "BasisIntegral"; case IndependentObjectEnum : return "IndependentObject"; case InversionControlParametersEnum : return "InversionControlParameters"; Index: /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp (revision 14966) +++ /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp (revision 14967) @@ -116,4 +116,5 @@ else if (strcmp(name,"HydrologyEfficient")==0) return HydrologyEfficientEnum; else if (strcmp(name,"HydrologySedimentKmax")==0) return HydrologySedimentKmaxEnum; + else if (strcmp(name,"BasisIntegral")==0) return BasisIntegralEnum; else if (strcmp(name,"IndependentObject")==0) return IndependentObjectEnum; else if (strcmp(name,"InversionControlParameters")==0) return InversionControlParametersEnum; @@ -136,9 +137,9 @@ else if (strcmp(name,"InversionVelObs")==0) return InversionVelObsEnum; else if (strcmp(name,"InversionVxObs")==0) return InversionVxObsEnum; - else if (strcmp(name,"InversionVyObs")==0) return InversionVyObsEnum; else stage=2; } if(stage==2){ - if (strcmp(name,"InversionVzObs")==0) return InversionVzObsEnum; + if (strcmp(name,"InversionVyObs")==0) return InversionVyObsEnum; + else if (strcmp(name,"InversionVzObs")==0) return InversionVzObsEnum; else if (strcmp(name,"MaskElementonfloatingice")==0) return MaskElementonfloatingiceEnum; else if (strcmp(name,"MaskElementongroundedice")==0) return MaskElementongroundediceEnum; @@ -259,9 +260,9 @@ else if (strcmp(name,"TransientIsgroundingline")==0) return TransientIsgroundinglineEnum; else if (strcmp(name,"TransientIsprognostic")==0) return TransientIsprognosticEnum; - else if (strcmp(name,"TransientIsthermal")==0) return TransientIsthermalEnum; else stage=3; } if(stage==3){ - if (strcmp(name,"TransientIsgia")==0) return TransientIsgiaEnum; + if (strcmp(name,"TransientIsthermal")==0) return TransientIsthermalEnum; + else if (strcmp(name,"TransientIsgia")==0) return TransientIsgiaEnum; else if (strcmp(name,"TransientNumRequestedOutputs")==0) return TransientNumRequestedOutputsEnum; else if (strcmp(name,"TransientRequestedOutputs")==0) return TransientRequestedOutputsEnum; @@ -382,9 +383,9 @@ else if (strcmp(name,"TriaP1Input")==0) return TriaP1InputEnum; else if (strcmp(name,"Vertex")==0) return VertexEnum; - else if (strcmp(name,"Air")==0) return AirEnum; else stage=4; } if(stage==4){ - if (strcmp(name,"Ice")==0) return IceEnum; + if (strcmp(name,"Air")==0) return AirEnum; + else if (strcmp(name,"Ice")==0) return IceEnum; else if (strcmp(name,"Melange")==0) return MelangeEnum; else if (strcmp(name,"Water")==0) return WaterEnum; @@ -505,9 +506,9 @@ else if (strcmp(name,"MinVel")==0) return MinVelEnum; else if (strcmp(name,"MaxVel")==0) return MaxVelEnum; - else if (strcmp(name,"MinVx")==0) return MinVxEnum; else stage=5; } if(stage==5){ - if (strcmp(name,"MaxVx")==0) return MaxVxEnum; + if (strcmp(name,"MinVx")==0) return MinVxEnum; + else if (strcmp(name,"MaxVx")==0) return MaxVxEnum; else if (strcmp(name,"MaxAbsVx")==0) return MaxAbsVxEnum; else if (strcmp(name,"MinVy")==0) return MinVyEnum; Index: /issm/trunk-jpl/src/c/solvers/solver_hydro_nonlinear.cpp =================================================================== --- /issm/trunk-jpl/src/c/solvers/solver_hydro_nonlinear.cpp (revision 14966) +++ /issm/trunk-jpl/src/c/solvers/solver_hydro_nonlinear.cpp (revision 14967) @@ -71,5 +71,5 @@ InputUpdateFromConstantx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,converged,ConvergedEnum); InputUpdateFromSolutionx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,ug); - /*InputUpdateFromConstantx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,sediment_kmax,HydrologySedimentKmaxEnum);*/ + InputUpdateFromConstantx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,sediment_kmax,HydrologySedimentKmaxEnum); break; } Index: /issm/trunk-jpl/src/m/enum/BasisIntegralEnum.m =================================================================== --- /issm/trunk-jpl/src/m/enum/BasisIntegralEnum.m (revision 14967) +++ /issm/trunk-jpl/src/m/enum/BasisIntegralEnum.m (revision 14967) @@ -0,0 +1,11 @@ +function macro=BasisIntegralEnum() +%BASISINTEGRALENUM - Enum of BasisIntegral +% +% WARNING: DO NOT MODIFY THIS FILE +% this file has been automatically generated by src/c/shared/Enum/Synchronize.sh +% Please read src/c/shared/Enum/README for more information +% +% Usage: +% macro=BasisIntegralEnum() + +macro=StringToEnum('BasisIntegral'); Index: /issm/trunk-jpl/src/m/enum/EnumDefinitions.py =================================================================== --- /issm/trunk-jpl/src/m/enum/EnumDefinitions.py (revision 14966) +++ /issm/trunk-jpl/src/m/enum/EnumDefinitions.py (revision 14967) @@ -1395,4 +1395,18 @@ return StringToEnum('HydrologySedimentKmax')[0] +def BasisIntegralEnum(): + """ + BASISINTEGRALENUM - Enum of BasisIntegral + + WARNING: DO NOT MODIFY THIS FILE + this file has been automatically generated by src/c/shared/Enum/Synchronize.sh + Please read src/c/shared/Enum/README for more information + + Usage: + macro=BasisIntegralEnum() + """ + + return StringToEnum('BasisIntegral')[0] + def IndependentObjectEnum(): """ @@ -7665,4 +7679,4 @@ """ - return 546 - + return 547 + Index: /issm/trunk-jpl/src/m/enum/MaximumNumberOfEnums.m =================================================================== --- /issm/trunk-jpl/src/m/enum/MaximumNumberOfEnums.m (revision 14966) +++ /issm/trunk-jpl/src/m/enum/MaximumNumberOfEnums.m (revision 14967) @@ -9,3 +9,3 @@ % macro=MaximumNumberOfEnums() -macro=546; +macro=547;