Index: /issm/trunk-jpl/src/c/analyses/HydrologyDCEfficientAnalysis.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/HydrologyDCEfficientAnalysis.cpp (revision 18974) +++ /issm/trunk-jpl/src/c/analyses/HydrologyDCEfficientAnalysis.cpp (revision 18975) @@ -24,10 +24,11 @@ iomodel->Constant(&isefficientlayer,HydrologydcIsefficientlayerEnum); + /*If not return*/ if(!isefficientlayer) return; + + /*If yes, initialize a flip flop counter*/ iomodel->FetchData(&eplflip_lock,HydrologydcEplflipLockEnum); - parameters->AddObject(new IntParam(HydrologydcEplflipLockEnum,eplflip_lock)); - /*Nothing for now*/ }/*}}}*/ void HydrologyDCEfficientAnalysis::UpdateElements(Elements* elements,IoModel* iomodel,int analysis_counter,int analysis_type){/*{{{*/ @@ -35,5 +36,4 @@ bool isefficientlayer; int hydrology_model; - /*Now, do we really want DC?*/ @@ -61,7 +61,6 @@ iomodel->FetchDataToInput(elements,HydrologydcEplInitialThicknessEnum); iomodel->FetchDataToInput(elements,HydrologydcEplMaxThicknessEnum); - iomodel->FetchDataToInput(elements,HydrologydcSedimentTransmitivityEnum); iomodel->FetchDataToInput(elements,HydrologydcEplThicknessEnum); - if(iomodel->domaintype!=Domain2DhorizontalEnum){ + if(iomodel->domaintype!=Domain2DhorizontalEnum){ iomodel->FetchDataToInput(elements,MeshVertexonbaseEnum); iomodel->FetchDataToInput(elements,MeshVertexonsurfaceEnum); @@ -85,4 +84,5 @@ iomodel->DeleteData(2,MeshVertexonbaseEnum,MeshVertexonsurfaceEnum); }/*}}}*/ + void HydrologyDCEfficientAnalysis::CreateConstraints(Constraints* constraints,IoModel* iomodel){/*{{{*/ @@ -98,6 +98,6 @@ IoModelToConstraintsx(constraints,iomodel,HydrologydcSpceplHeadEnum,HydrologyDCEfficientAnalysisEnum,P1Enum); - -}/*}}}*/ +}/*}}}*/ + void HydrologyDCEfficientAnalysis::CreateLoads(Loads* loads, IoModel* iomodel){/*{{{*/ /*Nothing for now*/ @@ -131,14 +131,14 @@ }/*}}}*/ - ElementVector* HydrologyDCEfficientAnalysis::CreateDVector(Element* element){/*{{{*/ +ElementVector* HydrologyDCEfficientAnalysis::CreateDVector(Element* element){/*{{{*/ /*Default, return NULL*/ return NULL; }/*}}}*/ - ElementMatrix* HydrologyDCEfficientAnalysis::CreateJacobianMatrix(Element* element){/*{{{*/ +ElementMatrix* HydrologyDCEfficientAnalysis::CreateJacobianMatrix(Element* element){/*{{{*/ _error_("Not implemented"); }/*}}}*/ - ElementMatrix* HydrologyDCEfficientAnalysis::CreateKMatrix(Element* element){/*{{{*/ +ElementMatrix* HydrologyDCEfficientAnalysis::CreateKMatrix(Element* element){/*{{{*/ /*Intermediaries*/ @@ -244,6 +244,6 @@ if(domaintype!=Domain2DhorizontalEnum){basalelement->DeleteMaterials(); delete basalelement;}; return Ke; - -}/*}}}*/ +}/*}}}*/ + ElementVector* HydrologyDCEfficientAnalysis::CreatePVector(Element* element){/*{{{*/ @@ -348,10 +348,13 @@ return pe; }/*}}}*/ + void HydrologyDCEfficientAnalysis::GetSolutionFromInputs(Vector* solution,Element* element){/*{{{*/ element->GetSolutionFromInputsOneDof(solution,EplHeadEnum); }/*}}}*/ + void HydrologyDCEfficientAnalysis::GradientJ(Vector* gradient,Element* element,int control_type,int control_index){/*{{{*/ _error_("Not implemented yet"); }/*}}}*/ + void HydrologyDCEfficientAnalysis::InputUpdateFromSolution(IssmDouble* solution,Element* element){/*{{{*/ @@ -375,5 +378,4 @@ /*Fetch number of nodes for this finite element*/ int numnodes = basalelement->GetNumberOfNodes(); - /*Fetch dof list and allocate solution vector*/ @@ -401,4 +403,5 @@ if(domaintype!=Domain2DhorizontalEnum){basalelement->DeleteMaterials(); delete basalelement;}; } /*}}}*/ + void HydrologyDCEfficientAnalysis::UpdateConstraints(FemModel* femmodel){/*{{{*/ /*Default, do nothing*/ @@ -415,4 +418,5 @@ return rho_freshwater*g*epl_porosity*(water_compressibility+(epl_compressibility/epl_porosity)); }/*}}}*/ + IssmDouble HydrologyDCEfficientAnalysis::SedimentStoring(Element* element){/*{{{*/ IssmDouble rho_freshwater = element->GetMaterialParameter(MaterialsRhoFreshwaterEnum); @@ -535,5 +539,5 @@ _error_("not Implemented Yet"); } - + int numnodes = element->GetNumberOfNodes(); IssmDouble* thickness = xNew(numnodes); @@ -602,6 +606,6 @@ xDelete(bed); } -} -/*}}}*/ +}/*}}}*/ + void HydrologyDCEfficientAnalysis::GetB(IssmDouble* B,Element* element,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/ /*Compute B matrix. B=[B1 B2 B3] where Bi is of size 3*NDOF2. @@ -631,4 +635,5 @@ xDelete(dbasis); }/*}}}*/ + void HydrologyDCEfficientAnalysis::HydrologyEPLGetMask(Vector* vec_mask, int* eplzigzag_counter, Element* element){ @@ -655,5 +660,5 @@ /*Intermediaries*/ - int penalty_lock; + int eplflip_lock; int numnodes =basalelement->GetNumberOfNodes(); IssmDouble* epl_thickness =xNew(numnodes); @@ -669,5 +674,5 @@ active_element_input->GetInputValue(&active_element); - basalelement->parameters->FindParam(&penalty_lock,HydrologydcEplflipLockEnum); + basalelement->parameters->FindParam(&eplflip_lock,HydrologydcEplflipLockEnum); basalelement-> GetInputListOnVertices(&old_active[0],HydrologydcMaskEplactiveNodeEnum); @@ -691,8 +696,9 @@ else if(old_active[i]>0.){ vec_mask->SetValue(basalelement->nodes[i]->Sid(),1.,INS_VAL); - /*If epl thickness gets under 10-3 initial thickness, close the layer*/ + /*If epl thickness gets under colapse thickness, close the layer*/ if(epl_thickness[i]nodes[i]->Lid()] ++; - if(eplzigzag_counter[basalelement->nodes[i]->Lid()]nodes[i]->Lid()]SetValue(basalelement->nodes[i]->Sid(),0.,INS_VAL); epl_thickness[i]=init_thick; @@ -705,5 +711,4 @@ /*Increase of the efficient system is needed if the epl head reach the maximum value (sediment max value for now)*/ GetHydrologyDCInefficientHmax(&h_max,basalelement,basalelement->nodes[i]); - if(eplhead[i]>=h_max && active_element){ for(j=0;j* active_vec, Element* element){ /*Constants*/ @@ -750,5 +756,6 @@ IssmDouble flag = 0.; IssmDouble* active = xNew(numnodes); - + + /*Pass the activity mask from elements to nodes*/ basalelement->GetInputListOnVertices(&active[0],HydrologydcMaskEplactiveNodeEnum); bool active_element; @@ -758,4 +765,5 @@ for(int i=0;i0.){ for(int i=0;iAddObject(new IntParam(HydrologydcPenaltyLockEnum,penalty_lock)); parameters->AddObject(new IntParam(HydrologydcMaxIterEnum,hydro_maxiter)); - -}/*}}}*/ +}/*}}}*/ + void HydrologyDCInefficientAnalysis::UpdateElements(Elements* elements,IoModel* iomodel,int analysis_counter,int analysis_type){/*{{{*/ @@ -97,19 +97,7 @@ if(isefficientlayer){ iomodel->FetchDataToInput(elements,HydrologydcMaskEplactiveNodeEnum); - - /* for(int i=0;iSize();i++){ */ - /* Element* element=xDynamicCast(elements->GetObjectByOffset(i)); */ - /* Input* node_mask_input = element->GetInput(HydrologydcMaskEplactiveNodeEnum); _assert_(node_mask_input); */ - - /* if(node_mask_input->Max()>0.) { */ - /* element_active = true; */ - /* } */ - /* else{ */ - /* element_active = false; */ - /* } */ - /* element->AddInput(new BoolInput(HydrologydcMaskEplactiveEltEnum,element_active)); */ - /* } */ - } -}/*}}}*/ + } +}/*}}}*/ + void HydrologyDCInefficientAnalysis::CreateNodes(Nodes* nodes,IoModel* iomodel){/*{{{*/ @@ -125,4 +113,5 @@ iomodel->DeleteData(2,MeshVertexonbaseEnum,MeshVertexonsurfaceEnum); }/*}}}*/ + void HydrologyDCInefficientAnalysis::CreateConstraints(Constraints* constraints,IoModel* iomodel){/*{{{*/ @@ -134,4 +123,5 @@ IoModelToConstraintsx(constraints,iomodel,HydrologydcSpcsedimentHeadEnum,HydrologyDCInefficientAnalysisEnum,P1Enum); }/*}}}*/ + void HydrologyDCInefficientAnalysis::CreateLoads(Loads* loads, IoModel* iomodel){/*{{{*/ @@ -165,11 +155,14 @@ _error_("not implemented"); }/*}}}*/ + ElementVector* HydrologyDCInefficientAnalysis::CreateDVector(Element* element){/*{{{*/ /*Default, return NULL*/ return NULL; }/*}}}*/ + ElementMatrix* HydrologyDCInefficientAnalysis::CreateJacobianMatrix(Element* element){/*{{{*/ _error_("Not implemented"); }/*}}}*/ + ElementMatrix* HydrologyDCInefficientAnalysis::CreateKMatrix(Element* element){/*{{{*/ @@ -275,4 +268,5 @@ return Ke; }/*}}}*/ + ElementVector* HydrologyDCInefficientAnalysis::CreatePVector(Element* element){/*{{{*/ @@ -312,5 +306,4 @@ IssmDouble* basis = xNew(numnodes); - /*Retrieve all inputs and parameters*/ basalelement->GetVerticesCoordinates(&xyz_list); @@ -342,7 +335,5 @@ /*Loading term*/ water_input->GetInputValue(&water_load,gauss); - scalar = Jdet*gauss->weight*(water_load); - if(dt!=0.) scalar = scalar*dt; for(int i=0;i(dbasis); }/*}}}*/ + void HydrologyDCInefficientAnalysis::GetSolutionFromInputs(Vector* solution,Element* element){/*{{{*/ element->GetSolutionFromInputsOneDof(solution,SedimentHeadEnum); }/*}}}*/ + void HydrologyDCInefficientAnalysis::GradientJ(Vector* gradient,Element* element,int control_type,int control_index){/*{{{*/ _error_("Not implemented yet"); }/*}}}*/ + void HydrologyDCInefficientAnalysis::InputUpdateFromSolution(IssmDouble* solution,Element* element){/*{{{*/ @@ -489,4 +484,5 @@ if(domaintype!=Domain2DhorizontalEnum){basalelement->DeleteMaterials(); delete basalelement;}; }/*}}}*/ + void HydrologyDCInefficientAnalysis::UpdateConstraints(FemModel* femmodel){/*{{{*/ /*Default, do nothing*/ @@ -551,4 +547,5 @@ return h_max; }/*}}}*/ + void HydrologyDCInefficientAnalysis::GetHydrologyDCInefficientHmax(IssmDouble* ph_max,Element* element, Node* innode){/*{{{*/ Index: /issm/trunk-jpl/src/c/classes/FemModel.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/FemModel.cpp (revision 18974) +++ /issm/trunk-jpl/src/c/classes/FemModel.cpp (revision 18975) @@ -1979,9 +1979,4 @@ counter=sum_counter; if(VerboseSolution()) _printf0_(" Number of active nodes L2 Projection: "<< counter <<"\n"); - - /*+++++++this is done in the solution sequence+++++++++++++++*/ - /*Update dof indexings*/ - // this->UpdateConstraintsx(); - /*++++++++++++++++++++++++++*/ -} -/*}}}*/ +} +/*}}}*/ Index: /issm/trunk-jpl/src/c/solutionsequences/solutionsequence_hydro_nonlinear.cpp =================================================================== --- /issm/trunk-jpl/src/c/solutionsequences/solutionsequence_hydro_nonlinear.cpp (revision 18974) +++ /issm/trunk-jpl/src/c/solutionsequences/solutionsequence_hydro_nonlinear.cpp (revision 18975) @@ -127,12 +127,4 @@ /*}}}*//*End of the sediment penalization loop*/ - /*Update EPL mask*/ - - /*++++++This is probably useless+++++++*/ - /* if(isefficientlayer){ */ - /* inefanalysis->ElementizeEplMask(femmodel); */ - /* } */ - /*+++++++++++++*/ - sedconverged=false;