Index: /issm/trunk-jpl/src/c/Makefile.am =================================================================== --- /issm/trunk-jpl/src/c/Makefile.am (revision 16684) +++ /issm/trunk-jpl/src/c/Makefile.am (revision 16685) @@ -465,4 +465,6 @@ ./analyses/HydrologyShreveAnalysis.h\ ./analyses/HydrologyShreveAnalysis.cpp\ + ./analyses/L2ProjectionEPLAnalysis.h\ + ./analyses/L2ProjectionEPLAnalysis.cpp\ ./cores/hydrology_core.cpp\ ./solutionsequences/solutionsequence_hydro_nonlinear.cpp Index: /issm/trunk-jpl/src/c/analyses/EnumToAnalysis.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/EnumToAnalysis.cpp (revision 16684) +++ /issm/trunk-jpl/src/c/analyses/EnumToAnalysis.cpp (revision 16685) @@ -20,4 +20,5 @@ case BalancevelocityAnalysisEnum : return new BalancevelocityAnalysis(); case L2ProjectionBaseAnalysisEnum : return new L2ProjectionBaseAnalysis(); + case L2ProjectionEPLAnalysisEnum : return new L2ProjectionEPLAnalysis(); case DamageEvolutionAnalysisEnum : return new DamageEvolutionAnalysis(); case StressbalanceAnalysisEnum : return new StressbalanceAnalysis(); Index: /issm/trunk-jpl/src/c/analyses/L2ProjectionEPLAnalysis.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/L2ProjectionEPLAnalysis.cpp (revision 16685) +++ /issm/trunk-jpl/src/c/analyses/L2ProjectionEPLAnalysis.cpp (revision 16685) @@ -0,0 +1,55 @@ +#include "./L2ProjectionEPLAnalysis.h" +#include "../toolkits/toolkits.h" +#include "../classes/classes.h" +#include "../shared/shared.h" +#include "../modules/modules.h" + +/*Model processing*/ +int L2ProjectionEPLAnalysis::DofsPerNode(int** doflist,int meshtype,int approximation){/*{{{*/ + return 1; +}/*}}}*/ +void L2ProjectionEPLAnalysis::UpdateParameters(Parameters* parameters,IoModel* iomodel,int solution_enum,int analysis_enum){/*{{{*/ +}/*}}}*/ +void L2ProjectionEPLAnalysis::UpdateElements(Elements* elements,IoModel* iomodel,int analysis_counter,int analysis_type){/*{{{*/ + + /*Update elements: */ + int counter=0; + for(int i=0;inumberofelements;i++){ + if(iomodel->my_elements[i]){ + Element* element=(Element*)elements->GetObjectByOffset(counter); + element->Update(i,iomodel,analysis_counter,analysis_type,P1Enum); + counter++; + } + } + + iomodel->FetchDataToInput(elements,EplHeadEnum); + iomodel->FetchDataToInput(elements,MaskIceLevelsetEnum); + if(iomodel->meshtype==Mesh3DEnum){ + iomodel->FetchDataToInput(elements,MeshElementonbedEnum); + } +}/*}}}*/ +void L2ProjectionEPLAnalysis::CreateNodes(Nodes* nodes,IoModel* iomodel){/*{{{*/ + + if(iomodel->meshtype==Mesh3DEnum){ + iomodel->FetchData(1,MeshVertexonbedEnum); + } + else if(iomodel->meshtype==Mesh2DverticalEnum){ + iomodel->FetchData(1,MeshVertexonbedEnum); + } + ::CreateNodes(nodes,iomodel,L2ProjectionEPLAnalysisEnum,P1Enum); + iomodel->DeleteData(1,MeshVertexonbedEnum); +}/*}}}*/ +void L2ProjectionEPLAnalysis::CreateConstraints(Constraints* constraints,IoModel* iomodel){/*{{{*/ + + /*No constraints*/ +}/*}}}*/ +void L2ProjectionEPLAnalysis::CreateLoads(Loads* loads, IoModel* iomodel){/*{{{*/ + + /*No loads*/ +}/*}}}*/ + +/*Numerics*/ +void L2ProjectionEPLAnalysis::GetSolutionFromInputs(Vector* solution,Element* element){/*{{{*/ + _error_("not implemented yet"); +}/*}}}*/ + Index: /issm/trunk-jpl/src/c/analyses/L2ProjectionEPLAnalysis.h =================================================================== --- /issm/trunk-jpl/src/c/analyses/L2ProjectionEPLAnalysis.h (revision 16685) +++ /issm/trunk-jpl/src/c/analyses/L2ProjectionEPLAnalysis.h (revision 16685) @@ -0,0 +1,22 @@ +/*! \file L2ProjectionEPLAnalysis.h + * \brief: header file for generic external result object + */ + +#ifndef _L2ProjectionEPLAnalysis_ +#define _L2ProjectionEPLAnalysis_ + +/*Headers*/ +#include "./Analysis.h" + +class L2ProjectionEPLAnalysis: public Analysis{ + + public: + int DofsPerNode(int** doflist,int meshtype,int approximation); + void UpdateParameters(Parameters* parameters,IoModel* iomodel,int solution_enum,int analysis_enum); + void UpdateElements(Elements* elements,IoModel* iomodel,int analysis_counter,int analysis_type); + void CreateNodes(Nodes* nodes,IoModel* iomodel); + void CreateConstraints(Constraints* constraints,IoModel* iomodel); + void CreateLoads(Loads* loads, IoModel* iomodel); + void GetSolutionFromInputs(Vector* solution,Element* element); +}; +#endif Index: /issm/trunk-jpl/src/c/analyses/analyses.h =================================================================== --- /issm/trunk-jpl/src/c/analyses/analyses.h (revision 16684) +++ /issm/trunk-jpl/src/c/analyses/analyses.h (revision 16685) @@ -32,6 +32,6 @@ #include "./StressbalanceVerticalAnalysis.h" #include "./L2ProjectionBaseAnalysis.h" +#include "./L2ProjectionEPLAnalysis.h" #include "./ThermalAnalysis.h" - -#include "EnumToAnalysis.h" +#include "./EnumToAnalysis.h" #endif Index: /issm/trunk-jpl/src/c/classes/Elements/Element.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Element.h (revision 16684) +++ /issm/trunk-jpl/src/c/classes/Elements/Element.h (revision 16685) @@ -154,4 +154,5 @@ virtual void HydrologyEPLGetActive(Vector* active)=0; virtual void ComputeEPLThickness(void)=0; + virtual void UpdateConstraintsL2ProjectionEPL(void)=0; #endif Index: /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp (revision 16684) +++ /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp (revision 16685) @@ -481,4 +481,7 @@ case HydrologyDCEfficientAnalysisEnum: return CreateKMatrixHydrologyDCEfficient(); + break; + case L2ProjectionEPLAnalysisEnum: + return CreateEPLDomainMassMatrix(); break; #endif @@ -697,4 +700,7 @@ return CreatePVectorHydrologyDCEfficient(); break; + case L2ProjectionEPLAnalysisEnum: + return CreatePVectorL2ProjectionEPL(); + break; #endif default: @@ -2258,4 +2264,8 @@ case HydrologyDCEfficientAnalysisEnum: InputUpdateFromSolutionOneDofCollapsed(solution,EplHeadEnum); + break; + case L2ProjectionEPLAnalysisEnum: + this->parameters->FindParam(&inputenum,InputToL2ProjectEnum); + InputUpdateFromSolutionOneDofCollapsed(solution,inputenum); break; #endif @@ -10622,4 +10632,17 @@ } /*}}}*/ +/*FUNCTION Penta::CreateEPLDomainMassMatrix {{{*/ +ElementMatrix* Penta::CreateEPLDomainMassMatrix(void){ + + if (!IsOnBed()) return NULL; + + Tria* tria=(Tria*)SpawnTria(0); //lower face is 0, upper face is 1. + ElementMatrix* Ke=tria->CreateEPLDomainMassMatrix(); + delete tria->material; delete tria; + + /*clean up and return*/ + return Ke; +} +/*}}}*/ /*FUNCTION Penta::CreatePVectorHydrologyDCInefficient {{{*/ ElementVector* Penta::CreatePVectorHydrologyDCInefficient(void){ @@ -10644,4 +10667,18 @@ Tria* tria=(Tria*)SpawnTria(0); //lower face is 0, upper face is 1. ElementVector* pe=tria->CreatePVectorHydrologyDCEfficient(); + delete tria->material; delete tria; + + /*Clean up and return*/ + return pe; +} +/*}}}*/ +/*FUNCTION Penta::CreatePVectorL@ProjectionEPL {{{*/ +ElementVector* Penta::CreatePVectorL2ProjectionEPL(void){ + + if (!IsOnBed()) return NULL; + + /*Call Tria function*/ + Tria* tria=(Tria*)SpawnTria(0); //lower face is 0, upper face is 1. + ElementVector* pe=tria->CreatePVectorL2ProjectionEPL(); delete tria->material; delete tria; @@ -10775,11 +10812,13 @@ /*Compute first the effective pressure in the EPL*/ EPL_N=gravity*((rho_ice*ice_thickness[i])-(rho_water*(eplhead[i]-bed[i]))); - + if(EPL_N<0.0)EPL_N=0.0; /*Get then the gradient of EPL heads*/ EPLgrad = epl_slopeX[i]+epl_slopeY[i]; /*And proceed to the real thing*/ - thickness[i] = old_thickness[i]*(1-((rho_water*gravity*dt)/(rho_ice* latentheat))*epl_conductivity*pow(EPLgrad,2.0)-2.0*(A*dt/(pow(n,n)))*EPL_N); + thickness[i] = old_thickness[i]*(1+((rho_water*gravity*dt)/(rho_ice*latentheat))*epl_conductivity*pow(EPLgrad,2.0)-2.0*(A*dt/(pow(n,n)))*(pow(EPL_N,n))); thickness[i+numdof2d]=thickness[i]; + printf("N, %e - thick term2, %e \n",EPL_N,-2.0*(A*dt/(pow(n,n)))*(pow(EPL_N,n))); + printf("old thick, %e - thick , %e \n",old_thickness[i],thickness[i]); } } @@ -10866,4 +10905,21 @@ /*Free ressources:*/ xDelete(doflist); +} +/*}}}*/ +/*FUNCTION Penta::UpdateConstraintsL2ProjectionEPL{{{*/ +void Penta::UpdateConstraintsL2ProjectionEPL(void){ + + IssmDouble activeEpl[NUMVERTICES]; + + + if(!IsOnBed()){ + for(int i=0;iNumberofNodes();i++)this->nodes[i]->Deactivate(); + } + else{ + GetInputListOnVertices(&activeEpl[0],HydrologydcMaskEplactiveEnum); + for(int i=0;i<3;i++){ + if(!activeEpl[i])this->nodes[i]->Deactivate(); + } + } } /*}}}*/ Index: /issm/trunk-jpl/src/c/classes/Elements/Penta.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Penta.h (revision 16684) +++ /issm/trunk-jpl/src/c/classes/Elements/Penta.h (revision 16685) @@ -325,12 +325,15 @@ ElementMatrix* CreateKMatrixHydrologyDCInefficient(void); ElementMatrix* CreateKMatrixHydrologyDCEfficient(void); + ElementMatrix* CreateEPLDomainMassMatrix(void); ElementVector* CreatePVectorHydrologyDCInefficient(void); ElementVector* CreatePVectorHydrologyDCEfficient(void); - void GetHydrologyDCInefficientHmax(IssmDouble* ph_max, Node* innode); - void GetHydrologyTransfer(Vector* transfer); - void HydrologyEPLGetActive(Vector* active_vec); - void HydrologyEPLGetMask(Vector* vec_mask); - void InputUpdateFromSolutionHydrologyDCInefficient(IssmDouble* solution); - void ComputeEPLThickness(void); + ElementVector* CreatePVectorL2ProjectionEPL(void); + void GetHydrologyDCInefficientHmax(IssmDouble* ph_max, Node* innode); + void GetHydrologyTransfer(Vector* transfer); + void HydrologyEPLGetActive(Vector* active_vec); + void HydrologyEPLGetMask(Vector* vec_mask); + void InputUpdateFromSolutionHydrologyDCInefficient(IssmDouble* solution); + void ComputeEPLThickness(void); + void UpdateConstraintsL2ProjectionEPL(void); #endif Index: /issm/trunk-jpl/src/c/classes/Elements/Seg.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Seg.h (revision 16684) +++ /issm/trunk-jpl/src/c/classes/Elements/Seg.h (revision 16685) @@ -117,4 +117,5 @@ void HydrologyEPLGetMask(Vector* vec_mask){_error_("not implemented yet");}; void ComputeEPLThickness(void){_error_("not implemented yet");}; + void UpdateConstraintsL2ProjectionEPL(){_error_("not implemented");}; #endif void GetSolutionFromInputsOneDof(Vector* solution,int enum_type){_error_("not implemented yet");}; Index: /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp (revision 16684) +++ /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp (revision 16685) @@ -288,4 +288,7 @@ case HydrologyDCEfficientAnalysisEnum: return CreateKMatrixHydrologyDCEfficient(); + break; + case L2ProjectionEPLAnalysisEnum: + return CreateEPLDomainMassMatrix(); break; #endif @@ -513,4 +516,7 @@ return CreatePVectorHydrologyDCEfficient(); break; + case L2ProjectionEPLAnalysisEnum: + return CreatePVectorL2ProjectionEPL(); + break; #endif #ifdef _HAVE_BALANCED_ @@ -586,6 +592,4 @@ case BedSlopeXEnum: input2 = inputs->GetInput(BedEnum); _assert_(input2); break; case BedSlopeYEnum: input2 = inputs->GetInput(BedEnum); _assert_(input2); break; - case EplHeadSlopeXEnum: input2 = inputs->GetInput(EplHeadEnum); _assert_(input2); break; - case EplHeadSlopeYEnum: input2 = inputs->GetInput(EplHeadEnum); _assert_(input2); break; default: input = inputs->GetInput(input_enum); } @@ -602,6 +606,6 @@ if(input2) input2->GetInputDerivativeValue(&slopes[0],&xyz_list[0][0],gauss); switch(input_enum){ - case SurfaceSlopeXEnum: case BedSlopeXEnum: case EplHeadSlopeXEnum: value = slopes[0]; break; - case SurfaceSlopeYEnum: case BedSlopeYEnum: case EplHeadSlopeYEnum: value = slopes[1]; break; + case SurfaceSlopeXEnum: case BedSlopeXEnum: value = slopes[0]; break; + case SurfaceSlopeYEnum: case BedSlopeYEnum: value = slopes[1]; break; default: input->GetInputValue(&value,gauss); } @@ -1784,4 +1788,8 @@ InputUpdateFromSolutionOneDof(solution,EplHeadEnum); break; + case L2ProjectionEPLAnalysisEnum: + this->parameters->FindParam(&extrusioninput,InputToL2ProjectEnum); + InputUpdateFromSolutionOneDof(solution,extrusioninput); + break; #endif #ifdef _HAVE_DAMAGE_ @@ -2901,5 +2909,15 @@ } /*}}}*/ - +/*FUNCTION Tria::UpdateConstraintsL2ProjectionEPL{{{*/ +void Tria::UpdateConstraintsL2ProjectionEPL(void){ + + IssmDouble activeEpl[NUMVERTICES]; + + GetInputListOnVertices(&activeEpl[0],HydrologydcMaskEplactiveEnum); + for(int i=0;inodes[i]->Deactivate(); + } +} +/*}}}*/ #ifdef _HAVE_RESPONSES_ /*FUNCTION Tria::AverageOntoPartition {{{*/ @@ -6683,4 +6701,44 @@ } /*}}}*/ + +/*FUNCTION Tria::CreatEPLDomainMassMatrix {{{*/ +ElementMatrix* Tria::CreateEPLDomainMassMatrix(void){ + + /* Intermediaries */ + IssmDouble D,Jdet; + IssmDouble xyz_list[NUMVERTICES][3]; + + /*Fetch number of nodes and dof for this finite element*/ + int numnodes = this->NumberofNodes(); + + /*Initialize Element matrix and vectors*/ + ElementMatrix* Ke = new ElementMatrix(nodes,numnodes,this->parameters,NoneApproximationEnum); + IssmDouble* basis = xNew(numnodes); + + /*Retrieve all inputs and parameters*/ + GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES); + + /* Start looping on the number of gaussian points: */ + GaussTria* gauss=new GaussTria(2); + for(int ig=gauss->begin();igend();ig++){ + + gauss->GaussPoint(ig); + + GetNodalFunctions(basis,gauss); + GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss); + D=gauss->weight*Jdet; + + TripleMultiply(basis,1,numnodes,1, + &D,1,1,0, + basis,1,numnodes,0, + &Ke->values[0],1); + } + + /*Clean up and return*/ + delete gauss; + xDelete(basis); + return Ke; +} +/*}}}*/ /*FUNCTION Tria::CreatePVectorHydrologyShreve {{{*/ ElementVector* Tria::CreatePVectorHydrologyShreve(void){ @@ -6708,5 +6766,4 @@ Input* basal_melting_input=inputs->GetInput(BasalforcingsMeltingRateEnum); _assert_(basal_melting_input); Input* old_watercolumn_input=inputs->GetInput(WaterColumnOldEnum); _assert_(old_watercolumn_input); - /*Initialize basal_melting_correction_g to 0, do not forget!:*/ /* Start looping on the number of gaussian points: */ @@ -6864,4 +6921,56 @@ residual_input->GetInputValue(&residual,gauss); pe->values[iv]+=residual/connectivity; + } + + /*Clean up and return*/ + xDelete(basis); + delete gauss; + return pe; +} +/*}}}*/ +/*FUNCTION Tria::CreatePVectorL2ProjectionEPL {{{*/ +ElementVector* Tria::CreatePVectorL2ProjectionEPL(void){ + + /*Intermediaries */ + int i,input_enum; + IssmDouble Jdet,value; + IssmDouble xyz_list[NUMVERTICES][3]; + IssmDouble slopes[2]; + Input* input = NULL; + Input* input2 = NULL; + + /*Fetch number of nodes and dof for this finite element*/ + int numnodes = this->NumberofNodes(); + + /*Initialize Element vector*/ + ElementVector* pe = new ElementVector(nodes,numnodes,this->parameters); + IssmDouble* basis = xNew(numnodes); + + /*Retrieve all inputs and parameters*/ + GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES); + this->parameters->FindParam(&input_enum,InputToL2ProjectEnum); + switch(input_enum){ + case EplHeadSlopeXEnum: input2 = inputs->GetInput(SurfaceEnum); _assert_(input2); break; + case EplHeadSlopeYEnum: input2 = inputs->GetInput(SurfaceEnum); _assert_(input2); break; + default: input = inputs->GetInput(input_enum); + } + + /* Start looping on the number of gaussian points: */ + GaussTria* gauss=new GaussTria(2); + for(int ig=gauss->begin();igend();ig++){ + + gauss->GaussPoint(ig); + + GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss); + GetNodalFunctions(basis,gauss); + + if(input2) input2->GetInputDerivativeValue(&slopes[0],&xyz_list[0][0],gauss); + switch(input_enum){ + case EplHeadSlopeXEnum: value = slopes[0]; break; + case EplHeadSlopeYEnum: value = slopes[1]; break; + default: input->GetInputValue(&value,gauss); + } + + for(i=0;ivalues[i]+=Jdet*gauss->weight*value*basis[i]; } @@ -7127,4 +7236,5 @@ IssmDouble h_max; IssmDouble sedheadmin; + IssmDouble epl_thickness[numdof]; IssmDouble old_active[numdof]; IssmDouble sedhead[numdof]; @@ -7133,4 +7243,5 @@ GetInputListOnVertices(&old_active[0],HydrologydcMaskEplactiveEnum); + GetInputListOnVertices(&epl_thickness[0],HydrologydcEplThicknessEnum); GetInputListOnVertices(&sedhead[0],SedimentHeadEnum); GetInputListOnVertices(&eplhead[0],EplHeadEnum); @@ -7149,7 +7260,11 @@ /*If mask was alread one, keep one*/ else if(old_active[i]>0.){ - vec_mask->SetValue(nodes[i]->Sid(),1.,INS_VAL); - } - + if(epl_thickness[i]>0.0){ + vec_mask->SetValue(nodes[i]->Sid(),1.,INS_VAL); + } + else{ + vec_mask->SetValue(nodes[i]->Sid(),0.,INS_VAL); + } + } /*Increase of the efficient system is needed if the epl head reach the maximum value (sediment max value for now)*/ this->GetHydrologyDCInefficientHmax(&h_max,nodes[i]); @@ -7217,11 +7332,11 @@ /*Compute first the effective pressure in the EPL*/ EPL_N=gravity*((rho_ice*ice_thickness[i])-(rho_water*(eplhead[i]-bed[i]))); - + if(EPL_N<0.0)EPL_N=0.0; /*Get then the gradient of EPL heads*/ EPLgrad = epl_slopeX[i]+epl_slopeY[i]; /*And proceed to the real thing*/ - thickness[i] = old_thickness[i]*(1-((rho_water*gravity*dt)/(rho_ice*latentheat))*epl_conductivity*pow(EPLgrad,2.0)-2.0*(A*dt/(pow(n,n)))*EPL_N); - + thickness[i] = old_thickness[i]*(1+((rho_water*gravity*dt)/(rho_ice*latentheat))*epl_conductivity*pow(EPLgrad,2.0)-2.0*(A*dt/(pow(n,n)))*(pow(EPL_N,n))); + if(this->nodes[i]->id==1553)printf("term1 %e, term2 %e\n",+((rho_water*gravity*dt)/(rho_ice*latentheat))*epl_conductivity*pow(EPLgrad,2.0), -2.0*(A*dt/(pow(n,n)))*(pow(EPL_N,n))); } } Index: /issm/trunk-jpl/src/c/classes/Elements/Tria.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Tria.h (revision 16684) +++ /issm/trunk-jpl/src/c/classes/Elements/Tria.h (revision 16685) @@ -301,20 +301,23 @@ ElementMatrix* CreateKMatrixHydrologyDCInefficient(void); ElementMatrix* CreateKMatrixHydrologyDCEfficient(void); + ElementMatrix* CreateEPLDomainMassMatrix(void); ElementVector* CreatePVectorHydrologyShreve(void); ElementVector* CreatePVectorHydrologyDCInefficient(void); ElementVector* CreatePVectorHydrologyDCEfficient(void); - void CreateHydrologyWaterVelocityInput(void); - void InputUpdateFromSolutionHydrology(IssmDouble* solution); - void InputUpdateFromSolutionHydrologyShreve(IssmDouble* solution); - void InputUpdateFromSolutionHydrologyDC(IssmDouble* solution); - void InputUpdateFromSolutionHydrologyDCInefficient(IssmDouble* solution); - void InputUpdateFromSolutionHydrologyDCEfficient(IssmDouble* solution); - void GetHydrologyDCInefficientHmax(IssmDouble* ph_max, Node* innode); - void GetHydrologyTransfer(Vector* transfer); - void HydrologyEPLGetActive(Vector* active_vec); - void HydrologyEPLGetMask(Vector* vec_mask); - void ComputeEPLThickness(void); - bool AllActive(void); - bool AnyActive(void); + ElementVector* CreatePVectorL2ProjectionEPL(void); + void CreateHydrologyWaterVelocityInput(void); + void InputUpdateFromSolutionHydrology(IssmDouble* solution); + void InputUpdateFromSolutionHydrologyShreve(IssmDouble* solution); + void InputUpdateFromSolutionHydrologyDC(IssmDouble* solution); + void InputUpdateFromSolutionHydrologyDCInefficient(IssmDouble* solution); + void InputUpdateFromSolutionHydrologyDCEfficient(IssmDouble* solution); + void GetHydrologyDCInefficientHmax(IssmDouble* ph_max, Node* innode); + void GetHydrologyTransfer(Vector* transfer); + void HydrologyEPLGetActive(Vector* active_vec); + void HydrologyEPLGetMask(Vector* vec_mask); + void ComputeEPLThickness(void); + void UpdateConstraintsL2ProjectionEPL(void); + bool AllActive(void); + bool AnyActive(void); #endif Index: /issm/trunk-jpl/src/c/classes/FemModel.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/FemModel.cpp (revision 16684) +++ /issm/trunk-jpl/src/c/classes/FemModel.cpp (revision 16685) @@ -261,7 +261,7 @@ * analyses. For example: do a SurfaceSlopeX, SurfaceSlopeY, BedSlopeX and BedSlopeY analysis using the * Slope configuration.*/ - int found=-1; for(int i=0;iSize();i++){ + Element* element=dynamic_cast(elements->GetObjectByOffset(i)); + element->UpdateConstraintsL2ProjectionEPL(); + } + +} +/*}}}*/#endif Index: /issm/trunk-jpl/src/c/classes/FemModel.h =================================================================== --- /issm/trunk-jpl/src/c/classes/FemModel.h (revision 16684) +++ /issm/trunk-jpl/src/c/classes/FemModel.h (revision 16685) @@ -103,4 +103,5 @@ void HydrologyEPLupdateDomainx(void); void HydrologyEPLThicknessx(void); + void UpdateConstraintsL2ProjectionEPLx(void); #endif }; Index: /issm/trunk-jpl/src/c/cores/AnalysisConfiguration.cpp =================================================================== --- /issm/trunk-jpl/src/c/cores/AnalysisConfiguration.cpp (revision 16684) +++ /issm/trunk-jpl/src/c/cores/AnalysisConfiguration.cpp (revision 16685) @@ -56,5 +56,5 @@ case HydrologySolutionEnum: - numanalyses=4; + numanalyses=5; analyses=xNew(numanalyses); analyses[0]=HydrologyShreveAnalysisEnum; @@ -62,4 +62,5 @@ analyses[2]=HydrologyDCEfficientAnalysisEnum; analyses[3]=L2ProjectionBaseAnalysisEnum; + analyses[4]=L2ProjectionEPLAnalysisEnum; break; Index: /issm/trunk-jpl/src/c/solutionsequences/solutionsequence_hydro_nonlinear.cpp =================================================================== --- /issm/trunk-jpl/src/c/solutionsequences/solutionsequence_hydro_nonlinear.cpp (revision 16684) +++ /issm/trunk-jpl/src/c/solutionsequences/solutionsequence_hydro_nonlinear.cpp (revision 16685) @@ -46,4 +46,5 @@ femmodel->parameters->FindParam(&eps_hyd,HydrologydcRelTolEnum); femmodel->parameters->FindParam(&time,TimeEnum); + /*FIXME, hardcoded, put on an enum*/ hydro_maxiter=150; hydrocount=1; @@ -98,4 +99,6 @@ if(num_unstable_constraints==0) sedconverged = true; if (sedcount>=hydro_maxiter){ + /*Hacking to get the results of non converged runs*/ + // sedconverged = true; _error_(" maximum number of Sediment iterations (" << hydro_maxiter << ") exceeded"); } @@ -116,10 +119,10 @@ /*Start by retrieving the EPL head slopes*/ - if(VerboseSolution()) _printf0_("computing EPL Head slope...\n"); - femmodel->SetCurrentConfiguration(L2ProjectionBaseAnalysisEnum); - femmodel->parameters->SetParam(EplHeadSlopeXEnum,InputToL2ProjectEnum); - solutionsequence_linear(femmodel); - femmodel->parameters->SetParam(EplHeadSlopeYEnum,InputToL2ProjectEnum); - solutionsequence_linear(femmodel); + /* if(VerboseSolution()) _printf0_("computing EPL Head slope...\n"); */ + /* femmodel->SetCurrentConfiguration(L2ProjectionEPLAnalysisEnum); */ + /* femmodel->parameters->SetParam(EplHeadSlopeXEnum,InputToL2ProjectEnum); */ + /* solutionsequence_linear(femmodel); */ + /* femmodel->parameters->SetParam(EplHeadSlopeYEnum,InputToL2ProjectEnum); */ + /* solutionsequence_linear(femmodel); */ femmodel->SetCurrentConfiguration(HydrologyDCEfficientAnalysisEnum); @@ -134,5 +137,18 @@ eplconverged = false; for(;;){ + + /*Start by retrieving the EPL head slopes*/ + if(VerboseSolution()) _printf0_("computing EPL Head slope...\n"); + femmodel->SetCurrentConfiguration(L2ProjectionEPLAnalysisEnum); + femmodel->UpdateConstraintsL2ProjectionEPLx(); + femmodel->parameters->SetParam(EplHeadSlopeXEnum,InputToL2ProjectEnum); + solutionsequence_linear(femmodel); + femmodel->parameters->SetParam(EplHeadSlopeYEnum,InputToL2ProjectEnum); + solutionsequence_linear(femmodel); + femmodel->SetCurrentConfiguration(HydrologyDCEfficientAnalysisEnum); + femmodel->HydrologyEPLThicknessx(); + //updating mask after the computation of the epl thickness + femmodel->HydrologyEPLupdateDomainx(); femmodel->HydrologyTransferx(); SystemMatricesx(&Kff,&Kfs,&pf,&df,NULL,femmodel); @@ -155,4 +171,6 @@ if(num_unstable_constraints==0) eplconverged = true; if (eplcount>=hydro_maxiter){ + /*Hacking to get the results of non converged runs*/ + //eplconverged = true; _error_(" maximum number of EPL iterations (" << hydro_maxiter << ") exceeded"); } @@ -214,5 +232,7 @@ else _printf0_(setw(50) << left << " Convergence criterion:" << ndu_sed/nu_sed*100 << " %\n"); if (hydrocount>=hydro_maxiter){ - _error_(" maximum number for hydrological global iterations (" << hydro_maxiter << ") exceeded"); + /*Hacking to get the results of non converged runs*/ + //hydroconverged = true; + _error_(" maximum number for hydrological global iterations (" << hydro_maxiter << ") exceeded"); } }