Index: /issm/trunk-jpl/src/c/classes/Elements/Element.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Element.h (revision 15248) +++ /issm/trunk-jpl/src/c/classes/Elements/Element.h (revision 15249) @@ -41,5 +41,4 @@ virtual void CreatePVector(Vector* pf)=0; virtual void CreateJacobianMatrix(Matrix* Jff)=0; - virtual void BasisIntegral(Vector* basisg)=0; virtual void GetSolutionFromInputs(Vector* solution)=0; virtual int GetNodeIndex(Node* node)=0; @@ -136,4 +135,6 @@ virtual void GetHydrologyTransfer(Vector* transfer)=0; virtual void UpdateConstraints(void)=0; + virtual void HydrologyEPLGetMask(Vector* mask)=0; + virtual void HydrologyEPLGetActive(Vector* active)=0; #endif }; Index: /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp (revision 15248) +++ /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp (revision 15249) @@ -802,40 +802,4 @@ /*return output*/ return penta; -} -/*}}}*/ -/*FUNCTION Penta::BasisIntegral {{{*/ -void Penta::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; - GaussPenta* gauss=NULL; - - GetDofList(&doflist,NoneApproximationEnum,GsetEnum); - GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES); - - /* Start looping on the number of gaussian points: */ - gauss=new GaussPenta(2,2); - for(int ig=gauss->begin();igend();ig++){ - - gauss->GaussPoint(ig); - - GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss); - GetNodalFunctionsP1(&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); } /*}}}*/ @@ -2166,5 +2130,4 @@ name==InversionVyObsEnum || name==InversionVzObsEnum || - name==BasisIntegralEnum || name==TemperatureEnum || name==EnthalpyEnum || @@ -9298,3 +9261,14 @@ } /*}}}*/ + +/*FUNCTION Penta::HydrologyEPLGetActive {{{*/ +void Penta::HydrologyEPLGetActive(Vector* active_vec){ + _error_("Hydrological stuff not suported in Penta"); +} +/*}}}*/ +/*FUNCTION Penta::HydrologyEPLGetMask{{{*/ +void Penta::HydrologyEPLGetMask(Vector* vec_mask){ + _error_("Hydrological stuff not suported in Penta"); +} +/*}}}*/ #endif Index: /issm/trunk-jpl/src/c/classes/Elements/Penta.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Penta.h (revision 15248) +++ /issm/trunk-jpl/src/c/classes/Elements/Penta.h (revision 15249) @@ -263,5 +263,4 @@ void InputUpdateFromSolutionDiagnosticVert( IssmDouble* solutiong); void InputUpdateFromSolutionDiagnosticStokes( IssmDouble* solutiong); - void BasisIntegral(Vector* gbasis); void GetSolutionFromInputsDiagnosticHoriz(Vector* solutiong); void GetSolutionFromInputsDiagnosticHutter(Vector* solutiong); @@ -307,4 +306,6 @@ void GetHydrologyDCInefficientHmax(IssmDouble* ph_max, Node* innode); void GetHydrologyTransfer(Vector* transfer); + void HydrologyEPLGetActive(Vector* active_vec); + void HydrologyEPLGetMask(Vector* vec_mask); #endif #ifdef _HAVE_THERMAL_ Index: /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp (revision 15248) +++ /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp (revision 15249) @@ -26,5 +26,5 @@ /*FUNCTION Tria::Tria(){{{*/ Tria::Tria(){ - + int i; @@ -1099,40 +1099,4 @@ _assert_(x2*y3 - y2*x3 + x1*y2 - y1*x2 + x3*y1 - y3*x1>0); return (x2*y3 - y2*x3 + x1*y2 - y1*x2 + x3*y1 - y3*x1)/2; -} -/*}}}*/ -/*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); - GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES); - - /* 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); } /*}}}*/ @@ -1818,5 +1782,5 @@ break; case HydrologyDCEfficientAnalysisEnum: - InputUpdateFromSolutionHydrologyDCEfficient(solution); + InputUpdateFromSolutionOneDof(solution,EplHeadEnum); break; #endif @@ -1849,8 +1813,7 @@ void Tria::InputUpdateFromSolutionOneDof(IssmDouble* solution,int enum_type){ - const int numdof = NDOF1*NUMVERTICES; - - int* doflist=NULL; - IssmDouble values[numdof]; + const int numdof = NDOF1*NUMVERTICES; + int* doflist = NULL; + IssmDouble values[numdof]; /*Get dof list: */ @@ -1981,4 +1944,16 @@ for(int i=0;i(values[i])) _error_("NaN found in solution vector"); + } + /*Add input to the element: */ + this->inputs->AddInput(new TriaP1Input(name,values)); + + /*Free ressources:*/ + xDelete(doflist); + return; + + case NodeSIdEnum: + for(int i=0;iSid()]; if(xIsNan(values[i])) _error_("NaN found in solution vector"); } @@ -2139,5 +2114,4 @@ name==HydrologydcMaskEplactiveEnum || name==WaterTransferEnum || - name==BasisIntegralEnum || name==QmuVxEnum || name==QmuVyEnum || @@ -2828,13 +2802,4 @@ /*Skip if water element*/ if(IsOnWater()) return; - - /*Just branch to the correct element stiffness matrix generator, according to the type of analysis we are carrying out: */ - switch(analysis_type){ - #ifdef _HAVE_HYDROLOGY_ - case HydrologyDCEfficientAnalysisEnum: - this->HydrologyUpdateEplDomain(); - break; - #endif - } } @@ -5829,15 +5794,27 @@ /*FUNCTION Tria::AllActive{{{*/ bool Tria::AllActive(void){ - + /*Intermediaries*/ const int numnodes = NUMVERTICES; - + for(int i=0;inodes[i]->IsActive()) return false; } - + return true; } /*}}}*/ +/*FUNCTION Tria::AnyActive{{{*/ +bool Tria::AnyActive(void){ + + /*Intermediaries*/ + const int numnodes = NUMVERTICES; + + for(int i=0;inodes[i]->IsActive()) return true; + } + + return false; +}/*}}}*/ /*FUNCTION Tria::CreateHydrologyWaterVelocityInput {{{*/ void Tria::CreateHydrologyWaterVelocityInput(void){ @@ -6076,5 +6053,4 @@ /*Check that all nodes are active, else return empty matrix*/ if(!this->AllActive()) return NULL; - /*Initialize Element matrix*/ ElementMatrix* Ke=new ElementMatrix(nodes,NUMVERTICES,this->parameters,NoneApproximationEnum); @@ -6500,28 +6476,4 @@ } /*}}}*/ -/*FUNCTION Tria::InputUpdateFromSolutionHydrologyDCEfficient{{{*/ -void Tria::InputUpdateFromSolutionHydrologyDCEfficient(IssmDouble* solution){ - - /*Intermediaries*/ - const int numdof = NDOF1 *NUMVERTICES; - int *doflist = NULL; - IssmDouble values[numdof]; - - /*Get dof list: */ - GetDofList(&doflist,NoneApproximationEnum,GsetEnum); - - /*Use the dof list to index into the solution vector: */ - for(int i=0;i(values[i])) _error_("NaN found in solution vector"); - } - - /*Add input to the element: */ - this->inputs->AddInput(new TriaP1Input(EplHeadEnum,values)); - - /*Free ressources:*/ - xDelete(doflist); -} -/*}}}*/ /*FUNCTION Tria::GetHydrologyDCInefficientHmax{{{*/ void Tria::GetHydrologyDCInefficientHmax(IssmDouble* ph_max, Node* innode){ @@ -6570,14 +6522,14 @@ IssmDouble sed_trans,sed_thick; IssmDouble leakage,h_max; + IssmDouble wh_trans; IssmDouble activeEpl[numdof]; - IssmDouble wh_trans[numdof]={0.0}; IssmDouble storing[numdof]; IssmDouble epl_head[numdof],sed_head[numdof]; GetDofList(&doflist,NoneApproximationEnum,GsetEnum); - + /*Get the flag to know if the efficient layer is present*/ - this->parameters->FindParam(&isefficientlayer,HydrologydcIsefficientlayerEnum); - + this->parameters->FindParam(&isefficientlayer,HydrologydcIsefficientlayerEnum); + if(isefficientlayer){ /*Also get the flag to the transfer method*/ @@ -6592,5 +6544,5 @@ GetInputListOnVertices(&activeEpl[0],HydrologydcMaskEplactiveEnum); - GetInputListOnVertices(&sed_head[0],SedimentHeadEnum); + GetInputListOnVertices(&sed_head[0],SedimentHeadEnum); GetInputListOnVertices(&epl_head[0],EplHeadEnum); @@ -6605,18 +6557,19 @@ if(sed_head[i]>epl_head[i]){ storing[i]=matpar->GetSedimentStoring(); - wh_trans[i]=sed_trans*storing[i]*(epl_head[i]-sed_head[i])/(leakage*sed_thick); + wh_trans=sed_trans*storing[i]*(epl_head[i]-sed_head[i])/(leakage*sed_thick); } else{ this->GetHydrologyDCInefficientHmax(&h_max,nodes[i]); storing[i]=matpar->GetEplStoring(); - wh_trans[i]=sed_trans*storing[i]*(epl_head[i]-sed_head[i])/(leakage*sed_thick); + wh_trans=sed_trans*storing[i]*(epl_head[i]-sed_head[i])/(leakage*sed_thick); if(sed_head[i]>h_max){ - wh_trans[i]=0.0; + wh_trans=0.0; } - if((sed_head[i]+wh_trans[i])>h_max){ - wh_trans[i]=h_max-sed_head[i]; + if((sed_head[i]+wh_trans)>h_max){ + wh_trans=h_max-sed_head[i]; } } - transfer->SetValue(doflist[i],wh_trans[i],INS_VAL); + /*Assign output pointer*/ + transfer->SetValue(doflist[i],wh_trans,INS_VAL); } break; @@ -6625,29 +6578,73 @@ } } - /*Assign output pointer*/ - /*transfer->SetValues(numdof,doflist,&wh_trans[0],INS_VAL);*/ -} - -/*}}}*/ -/*FUNCTION Tria::HydrologyUpdateEplDomain{{{*/ -void Tria::HydrologyUpdateEplDomain(void){ +} +/*}}}*/ +/*FUNCTION Tria::HydrologyEPLGetActive {{{*/ +void Tria::HydrologyEPLGetActive(Vector* active_vec){ + + /*Constants*/ + const int numnodes =NUMVERTICES; + int flag; + IssmDouble active[numnodes]; + + GetInputListOnVertices(&active[0],HydrologydcMaskEplactiveEnum); + + for(int i=0;i0.){ + for(int i=0;iActivate(); + } + + else if(!this->AnyActive()){ + for(int i=0;iDeactivate(); + + } + else{ + /*Do not do anything: at least one node is active for this element but this element is not solved for*/ + } + +} +/*}}}*/ +/*FUNCTION Tria::HydrologyEPLGetMask{{{*/ +void Tria::HydrologyEPLGetMask(Vector* vec_mask){ /*Intermediaries*/ - const int numdof = NDOF1 *NUMVERTICES; - bool isefficientlayer; - IssmDouble activeEpl[numdof]; - - /*Get parameter and Inout list*/ - this->parameters->FindParam(&isefficientlayer,HydrologydcIsefficientlayerEnum); - GetInputListOnVertices(&activeEpl[0],HydrologydcMaskEplactiveEnum); - - if(isefficientlayer){ - for(int i=0;iActivate(); - else if (activeEpl[i]==0.0) - nodes[i]->Deactivate(); - else - _error_("activation value "<< activeEpl[i] <<" not supported"); + int i,j; + const int numdof = NDOF1 *NUMVERTICES; + IssmDouble h_max; + IssmDouble sedheadmin; + IssmDouble old_active[numdof]; + IssmDouble sedhead[numdof]; + IssmDouble eplhead[numdof]; + + + GetInputListOnVertices(&old_active[0],HydrologydcMaskEplactiveEnum); + GetInputListOnVertices(&sedhead[0],SedimentHeadEnum); + GetInputListOnVertices(&eplhead[0],EplHeadEnum); + + sedheadmin=sedhead[0]; + for(i=1;i0.){ + vec_mask->SetValue(nodes[i]->Sid(),1.,INS_VAL); + } + + this->GetHydrologyDCInefficientHmax(&h_max,nodes[i]); + /*Increase the size of the efficient system if needed*/ + /*Increase is needed if the epl head reach the maximum value (sediment value for now)*/ + if(eplhead[i]>=h_max){ + vec_mask->SetValue(nodes[i]->Sid(),1.,INS_VAL); + for(j=0;jSetValue(nodes[j]->Sid(),1.,INS_VAL); + break; + } + } } } Index: /issm/trunk-jpl/src/c/classes/Elements/Tria.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Tria.h (revision 15248) +++ /issm/trunk-jpl/src/c/classes/Elements/Tria.h (revision 15249) @@ -89,5 +89,5 @@ int Sid(); bool IsOnBed(); - bool IsFloating(); + bool IsFloating(); bool IsNodeOnShelf(); bool IsNodeOnShelfFromFlags(IssmDouble* flags); @@ -196,5 +196,4 @@ ElementVector* CreatePVectorSlope(void); IssmDouble GetArea(void); - void BasisIntegral(Vector* gbasis); int GetElementType(void); void GetDofList(int** pdoflist,int approximation_enum,int setenum); @@ -250,5 +249,4 @@ void GetSolutionFromInputsHydrologyDCEfficient(Vector* solution); void CreateHydrologyWaterVelocityInput(void); - void HydrologyUpdateEplDomain(void); void UpdateConstraints(void); void InputUpdateFromSolutionHydrology(IssmDouble* solution); @@ -259,5 +257,8 @@ void GetHydrologyDCInefficientHmax(IssmDouble* ph_max, Node* innode); void GetHydrologyTransfer(Vector* transfer); + void HydrologyEPLGetActive(Vector* active_vec); + void HydrologyEPLGetMask(Vector* vec_mask); bool AllActive(void); + bool AnyActive(void); #endif #ifdef _HAVE_BALANCED_ Index: /issm/trunk-jpl/src/c/classes/FemModel.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/FemModel.cpp (revision 15248) +++ /issm/trunk-jpl/src/c/classes/FemModel.cpp (revision 15249) @@ -1581,20 +1581,50 @@ #endif -void FemModel::BasisIntegralsx(void){ /*{{{*/ - - Vector* basisg=NULL; - - /*Vector allocation*/ - basisg=new Vector(nodes->NumberOfNodes()); - +void FemModel::HydrologyEPLupdateDomainx(void){ /*{{{*/ + + Vector *mask = NULL; + IssmDouble *serial_mask = NULL; + Vector *active = NULL; + IssmDouble *serial_active = NULL; + + /*Step 1: update maks, the mask might be extended by residual and/or using downstream sediment head*/ + mask=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); + element->HydrologyEPLGetMask(mask); + } + + /*Assemble and serialize*/ + mask->Assemble(); + serial_mask=mask->ToMPISerial(); + delete mask; + + /*Update Mask*/ + InputUpdateFromVectorx(elements,nodes,vertices,loads,materials,parameters,serial_mask,HydrologydcMaskEplactiveEnum,NodeSIdEnum); + xDelete(serial_mask); + + /*Step 2: update node activity. If one element is connected to mask=1, all nodes are active*/ + active=new Vector(nodes->NumberOfNodes()); + for (int i=0;iSize();i++){ + Element* element=dynamic_cast(elements->GetObjectByOffset(i)); + element->HydrologyEPLGetActive(active); + } + + /*Assemble and serialize*/ + active->Assemble(); + serial_active=active->ToMPISerial(); + delete active; + + /*Update node activation accordingly*/ + for (int i=0;iSize();i++){ + Node* node=dynamic_cast(nodes->GetObjectByOffset(i)); + if(serial_active[node->Sid()]==1.){ + node->Activate(); + } + else{ + node->Deactivate(); + } + } + xDelete(serial_active); } Index: /issm/trunk-jpl/src/c/classes/FemModel.h =================================================================== --- /issm/trunk-jpl/src/c/classes/FemModel.h (revision 15248) +++ /issm/trunk-jpl/src/c/classes/FemModel.h (revision 15249) @@ -99,6 +99,7 @@ void UpdateConstraintsx(void); int UpdateVertexPositionsx(void); - void BasisIntegralsx(void); + void ParEplMask(void); void HydrologyTransferx(void); + void HydrologyEPLupdateDomainx(void); }; Index: /issm/trunk-jpl/src/c/solutionsequences/solutionsequence_hydro_nonlinear.cpp =================================================================== --- /issm/trunk-jpl/src/c/solutionsequences/solutionsequence_hydro_nonlinear.cpp (revision 15248) +++ /issm/trunk-jpl/src/c/solutionsequences/solutionsequence_hydro_nonlinear.cpp (revision 15249) @@ -41,5 +41,4 @@ femmodel->parameters->FindParam(&eps_hyd,HydrologydcRelTolEnum); femmodel->parameters->FindParam(&time,TimeEnum); - femmodel->BasisIntegralsx(); hydro_maxiter=100; hydrocount=1; @@ -123,5 +122,5 @@ delete uf; Solverx(&uf, Kff, pf,old_uf, df, femmodel->parameters); - delete old_uf; old_uf=uf->Duplicate(); + delete old_uf; old_uf=uf->Duplicate(); if(eplcount>1) delete ug_epl; delete Kff;delete pf; @@ -129,4 +128,5 @@ Mergesolutionfromftogx(&ug_epl,uf,ys,femmodel->nodes,femmodel->parameters); delete ys; InputUpdateFromSolutionx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,ug_epl); + femmodel->HydrologyEPLupdateDomainx(); ConstraintsStatex(&constraints_converged,&num_unstable_constraints, femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters); @@ -144,4 +144,5 @@ InputUpdateFromConstantx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,sediment_kmax,MeltingOffsetEnum); InputUpdateFromSolutionx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,ug_epl); + femmodel->HydrologyEPLupdateDomainx(); break; }