Changeset 15249
- Timestamp:
- 06/13/13 08:36:55 (12 years ago)
- Location:
- issm/trunk-jpl/src/c
- Files:
-
- 8 edited
Legend:
- Unmodified
- Added
- Removed
-
issm/trunk-jpl/src/c/classes/Elements/Element.h
r15223 r15249 41 41 virtual void CreatePVector(Vector<IssmDouble>* pf)=0; 42 42 virtual void CreateJacobianMatrix(Matrix<IssmDouble>* Jff)=0; 43 virtual void BasisIntegral(Vector<IssmDouble>* basisg)=0;44 43 virtual void GetSolutionFromInputs(Vector<IssmDouble>* solution)=0; 45 44 virtual int GetNodeIndex(Node* node)=0; … … 136 135 virtual void GetHydrologyTransfer(Vector<IssmDouble>* transfer)=0; 137 136 virtual void UpdateConstraints(void)=0; 137 virtual void HydrologyEPLGetMask(Vector<IssmDouble>* mask)=0; 138 virtual void HydrologyEPLGetActive(Vector<IssmDouble>* active)=0; 138 139 #endif 139 140 }; -
issm/trunk-jpl/src/c/classes/Elements/Penta.cpp
r15223 r15249 802 802 /*return output*/ 803 803 return penta; 804 }805 /*}}}*/806 /*FUNCTION Penta::BasisIntegral {{{*/807 void Penta::BasisIntegral(Vector<IssmDouble>* basisg){808 809 /*Constants*/810 const int numdof=NDOF1*NUMVERTICES;811 812 /*Intermediaries */813 IssmDouble Jdet;814 IssmDouble xyz_list[NUMVERTICES][3];815 IssmDouble basis[numdof];816 IssmDouble basisint[numdof]={0.};817 int *doflist=NULL;818 GaussPenta* gauss=NULL;819 820 GetDofList(&doflist,NoneApproximationEnum,GsetEnum);821 GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);822 823 /* Start looping on the number of gaussian points: */824 gauss=new GaussPenta(2,2);825 for(int ig=gauss->begin();ig<gauss->end();ig++){826 827 gauss->GaussPoint(ig);828 829 GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss);830 GetNodalFunctionsP1(&basis[0], gauss);831 832 for(int i=0;i<numdof;i++) basisint[i]+=Jdet*gauss->weight*basis[i];833 }834 835 basisg->SetValues(numdof,doflist,&basisint[0],ADD_VAL);836 837 /*Clean up and return*/838 delete gauss;839 xDelete<int>(doflist);840 804 } 841 805 /*}}}*/ … … 2166 2130 name==InversionVyObsEnum || 2167 2131 name==InversionVzObsEnum || 2168 name==BasisIntegralEnum ||2169 2132 name==TemperatureEnum || 2170 2133 name==EnthalpyEnum || … … 9298 9261 } 9299 9262 /*}}}*/ 9263 9264 /*FUNCTION Penta::HydrologyEPLGetActive {{{*/ 9265 void Penta::HydrologyEPLGetActive(Vector<IssmDouble>* active_vec){ 9266 _error_("Hydrological stuff not suported in Penta"); 9267 } 9268 /*}}}*/ 9269 /*FUNCTION Penta::HydrologyEPLGetMask{{{*/ 9270 void Penta::HydrologyEPLGetMask(Vector<IssmDouble>* vec_mask){ 9271 _error_("Hydrological stuff not suported in Penta"); 9272 } 9273 /*}}}*/ 9300 9274 #endif -
issm/trunk-jpl/src/c/classes/Elements/Penta.h
r15223 r15249 263 263 void InputUpdateFromSolutionDiagnosticVert( IssmDouble* solutiong); 264 264 void InputUpdateFromSolutionDiagnosticStokes( IssmDouble* solutiong); 265 void BasisIntegral(Vector<IssmDouble>* gbasis);266 265 void GetSolutionFromInputsDiagnosticHoriz(Vector<IssmDouble>* solutiong); 267 266 void GetSolutionFromInputsDiagnosticHutter(Vector<IssmDouble>* solutiong); … … 307 306 void GetHydrologyDCInefficientHmax(IssmDouble* ph_max, Node* innode); 308 307 void GetHydrologyTransfer(Vector<IssmDouble>* transfer); 308 void HydrologyEPLGetActive(Vector<IssmDouble>* active_vec); 309 void HydrologyEPLGetMask(Vector<IssmDouble>* vec_mask); 309 310 #endif 310 311 #ifdef _HAVE_THERMAL_ -
issm/trunk-jpl/src/c/classes/Elements/Tria.cpp
r15234 r15249 26 26 /*FUNCTION Tria::Tria(){{{*/ 27 27 Tria::Tria(){ 28 28 29 29 int i; 30 30 … … 1099 1099 _assert_(x2*y3 - y2*x3 + x1*y2 - y1*x2 + x3*y1 - y3*x1>0); 1100 1100 return (x2*y3 - y2*x3 + x1*y2 - y1*x2 + x3*y1 - y3*x1)/2; 1101 }1102 /*}}}*/1103 /*FUNCTION Tria::BasisIntegral {{{*/1104 void Tria::BasisIntegral(Vector<IssmDouble>* basisg){1105 1106 /*Constants*/1107 const int numdof=NDOF1*NUMVERTICES;1108 1109 /*Intermediaries */1110 IssmDouble Jdet;1111 IssmDouble xyz_list[NUMVERTICES][3];1112 IssmDouble basis[numdof];1113 IssmDouble basisint[numdof]={0.};1114 int *doflist=NULL;1115 GaussTria* gauss=NULL;1116 1117 GetDofList(&doflist,NoneApproximationEnum,GsetEnum);1118 GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);1119 1120 /* Start looping on the number of gaussian points: */1121 gauss=new GaussTria(2);1122 for(int ig=gauss->begin();ig<gauss->end();ig++){1123 1124 gauss->GaussPoint(ig);1125 1126 GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss);1127 GetNodalFunctions(&basis[0], gauss);1128 1129 for(int i=0;i<numdof;i++) basisint[i]+=Jdet*gauss->weight*basis[i];1130 }1131 1132 basisg->SetValues(numdof,doflist,&basisint[0],ADD_VAL);1133 1134 /*Clean up and return*/1135 delete gauss;1136 xDelete<int>(doflist);1137 1101 } 1138 1102 /*}}}*/ … … 1818 1782 break; 1819 1783 case HydrologyDCEfficientAnalysisEnum: 1820 InputUpdateFromSolution HydrologyDCEfficient(solution);1784 InputUpdateFromSolutionOneDof(solution,EplHeadEnum); 1821 1785 break; 1822 1786 #endif … … 1849 1813 void Tria::InputUpdateFromSolutionOneDof(IssmDouble* solution,int enum_type){ 1850 1814 1851 const int numdof = NDOF1*NUMVERTICES; 1852 1853 int* doflist=NULL; 1854 IssmDouble values[numdof]; 1815 const int numdof = NDOF1*NUMVERTICES; 1816 int* doflist = NULL; 1817 IssmDouble values[numdof]; 1855 1818 1856 1819 /*Get dof list: */ … … 1981 1944 for(int i=0;i<numdof;i++){ 1982 1945 values[i]=vector[doflist[i]]; 1946 if(xIsNan<IssmDouble>(values[i])) _error_("NaN found in solution vector"); 1947 } 1948 /*Add input to the element: */ 1949 this->inputs->AddInput(new TriaP1Input(name,values)); 1950 1951 /*Free ressources:*/ 1952 xDelete<int>(doflist); 1953 return; 1954 1955 case NodeSIdEnum: 1956 for(int i=0;i<NUMVERTICES;i++){ 1957 values[i]=vector[nodes[i]->Sid()]; 1983 1958 if(xIsNan<IssmDouble>(values[i])) _error_("NaN found in solution vector"); 1984 1959 } … … 2139 2114 name==HydrologydcMaskEplactiveEnum || 2140 2115 name==WaterTransferEnum || 2141 name==BasisIntegralEnum ||2142 2116 name==QmuVxEnum || 2143 2117 name==QmuVyEnum || … … 2828 2802 /*Skip if water element*/ 2829 2803 if(IsOnWater()) return; 2830 2831 /*Just branch to the correct element stiffness matrix generator, according to the type of analysis we are carrying out: */2832 switch(analysis_type){2833 #ifdef _HAVE_HYDROLOGY_2834 case HydrologyDCEfficientAnalysisEnum:2835 this->HydrologyUpdateEplDomain();2836 break;2837 #endif2838 }2839 2804 2840 2805 } … … 5829 5794 /*FUNCTION Tria::AllActive{{{*/ 5830 5795 bool Tria::AllActive(void){ 5831 5796 5832 5797 /*Intermediaries*/ 5833 5798 const int numnodes = NUMVERTICES; 5834 5799 5835 5800 for(int i=0;i<numnodes;i++){ 5836 5801 if(!this->nodes[i]->IsActive()) return false; 5837 5802 } 5838 5803 5839 5804 return true; 5840 5805 } 5841 5806 /*}}}*/ 5807 /*FUNCTION Tria::AnyActive{{{*/ 5808 bool Tria::AnyActive(void){ 5809 5810 /*Intermediaries*/ 5811 const int numnodes = NUMVERTICES; 5812 5813 for(int i=0;i<numnodes;i++){ 5814 if(this->nodes[i]->IsActive()) return true; 5815 } 5816 5817 return false; 5818 }/*}}}*/ 5842 5819 /*FUNCTION Tria::CreateHydrologyWaterVelocityInput {{{*/ 5843 5820 void Tria::CreateHydrologyWaterVelocityInput(void){ … … 6076 6053 /*Check that all nodes are active, else return empty matrix*/ 6077 6054 if(!this->AllActive()) return NULL; 6078 6079 6055 /*Initialize Element matrix*/ 6080 6056 ElementMatrix* Ke=new ElementMatrix(nodes,NUMVERTICES,this->parameters,NoneApproximationEnum); … … 6500 6476 } 6501 6477 /*}}}*/ 6502 /*FUNCTION Tria::InputUpdateFromSolutionHydrologyDCEfficient{{{*/6503 void Tria::InputUpdateFromSolutionHydrologyDCEfficient(IssmDouble* solution){6504 6505 /*Intermediaries*/6506 const int numdof = NDOF1 *NUMVERTICES;6507 int *doflist = NULL;6508 IssmDouble values[numdof];6509 6510 /*Get dof list: */6511 GetDofList(&doflist,NoneApproximationEnum,GsetEnum);6512 6513 /*Use the dof list to index into the solution vector: */6514 for(int i=0;i<numdof;i++){6515 values[i]=solution[doflist[i]];6516 if(xIsNan<IssmDouble>(values[i])) _error_("NaN found in solution vector");6517 }6518 6519 /*Add input to the element: */6520 this->inputs->AddInput(new TriaP1Input(EplHeadEnum,values));6521 6522 /*Free ressources:*/6523 xDelete<int>(doflist);6524 }6525 /*}}}*/6526 6478 /*FUNCTION Tria::GetHydrologyDCInefficientHmax{{{*/ 6527 6479 void Tria::GetHydrologyDCInefficientHmax(IssmDouble* ph_max, Node* innode){ … … 6570 6522 IssmDouble sed_trans,sed_thick; 6571 6523 IssmDouble leakage,h_max; 6524 IssmDouble wh_trans; 6572 6525 IssmDouble activeEpl[numdof]; 6573 IssmDouble wh_trans[numdof]={0.0};6574 6526 IssmDouble storing[numdof]; 6575 6527 IssmDouble epl_head[numdof],sed_head[numdof]; 6576 6528 6577 6529 GetDofList(&doflist,NoneApproximationEnum,GsetEnum); 6578 6530 6579 6531 /*Get the flag to know if the efficient layer is present*/ 6580 6581 6532 this->parameters->FindParam(&isefficientlayer,HydrologydcIsefficientlayerEnum); 6533 6582 6534 if(isefficientlayer){ 6583 6535 /*Also get the flag to the transfer method*/ … … 6592 6544 6593 6545 GetInputListOnVertices(&activeEpl[0],HydrologydcMaskEplactiveEnum); 6594 GetInputListOnVertices(&sed_head[0],SedimentHeadEnum); 6546 GetInputListOnVertices(&sed_head[0],SedimentHeadEnum); 6595 6547 GetInputListOnVertices(&epl_head[0],EplHeadEnum); 6596 6548 … … 6605 6557 if(sed_head[i]>epl_head[i]){ 6606 6558 storing[i]=matpar->GetSedimentStoring(); 6607 wh_trans [i]=sed_trans*storing[i]*(epl_head[i]-sed_head[i])/(leakage*sed_thick);6559 wh_trans=sed_trans*storing[i]*(epl_head[i]-sed_head[i])/(leakage*sed_thick); 6608 6560 } 6609 6561 else{ 6610 6562 this->GetHydrologyDCInefficientHmax(&h_max,nodes[i]); 6611 6563 storing[i]=matpar->GetEplStoring(); 6612 wh_trans [i]=sed_trans*storing[i]*(epl_head[i]-sed_head[i])/(leakage*sed_thick);6564 wh_trans=sed_trans*storing[i]*(epl_head[i]-sed_head[i])/(leakage*sed_thick); 6613 6565 if(sed_head[i]>h_max){ 6614 wh_trans [i]=0.0;6566 wh_trans=0.0; 6615 6567 } 6616 if((sed_head[i]+wh_trans [i])>h_max){6617 wh_trans [i]=h_max-sed_head[i];6568 if((sed_head[i]+wh_trans)>h_max){ 6569 wh_trans=h_max-sed_head[i]; 6618 6570 } 6619 6571 } 6620 transfer->SetValue(doflist[i],wh_trans[i],INS_VAL); 6572 /*Assign output pointer*/ 6573 transfer->SetValue(doflist[i],wh_trans,INS_VAL); 6621 6574 } 6622 6575 break; … … 6625 6578 } 6626 6579 } 6627 /*Assign output pointer*/ 6628 /*transfer->SetValues(numdof,doflist,&wh_trans[0],INS_VAL);*/ 6629 } 6630 6631 /*}}}*/ 6632 /*FUNCTION Tria::HydrologyUpdateEplDomain{{{*/ 6633 void Tria::HydrologyUpdateEplDomain(void){ 6580 } 6581 /*}}}*/ 6582 /*FUNCTION Tria::HydrologyEPLGetActive {{{*/ 6583 void Tria::HydrologyEPLGetActive(Vector<IssmDouble>* active_vec){ 6584 6585 /*Constants*/ 6586 const int numnodes =NUMVERTICES; 6587 int flag; 6588 IssmDouble active[numnodes]; 6589 6590 GetInputListOnVertices(&active[0],HydrologydcMaskEplactiveEnum); 6591 6592 for(int i=0;i<numnodes;i++) flag+=active[i]; 6593 6594 if(flag>0.){ 6595 for(int i=0;i<numnodes;i++) nodes[i]->Activate(); 6596 } 6597 6598 else if(!this->AnyActive()){ 6599 for(int i=0;i<numnodes;i++) nodes[i]->Deactivate(); 6600 6601 } 6602 else{ 6603 /*Do not do anything: at least one node is active for this element but this element is not solved for*/ 6604 } 6605 6606 } 6607 /*}}}*/ 6608 /*FUNCTION Tria::HydrologyEPLGetMask{{{*/ 6609 void Tria::HydrologyEPLGetMask(Vector<IssmDouble>* vec_mask){ 6634 6610 6635 6611 /*Intermediaries*/ 6636 const int numdof = NDOF1 *NUMVERTICES; 6637 bool isefficientlayer; 6638 IssmDouble activeEpl[numdof]; 6639 6640 /*Get parameter and Inout list*/ 6641 this->parameters->FindParam(&isefficientlayer,HydrologydcIsefficientlayerEnum); 6642 GetInputListOnVertices(&activeEpl[0],HydrologydcMaskEplactiveEnum); 6643 6644 if(isefficientlayer){ 6645 for(int i=0;i<numdof;i++){ 6646 if(activeEpl[i]==1.0) 6647 nodes[i]->Activate(); 6648 else if (activeEpl[i]==0.0) 6649 nodes[i]->Deactivate(); 6650 else 6651 _error_("activation value "<< activeEpl[i] <<" not supported"); 6612 int i,j; 6613 const int numdof = NDOF1 *NUMVERTICES; 6614 IssmDouble h_max; 6615 IssmDouble sedheadmin; 6616 IssmDouble old_active[numdof]; 6617 IssmDouble sedhead[numdof]; 6618 IssmDouble eplhead[numdof]; 6619 6620 6621 GetInputListOnVertices(&old_active[0],HydrologydcMaskEplactiveEnum); 6622 GetInputListOnVertices(&sedhead[0],SedimentHeadEnum); 6623 GetInputListOnVertices(&eplhead[0],EplHeadEnum); 6624 6625 sedheadmin=sedhead[0]; 6626 for(i=1;i<numdof;i++){ 6627 if(sedhead[i]<=sedheadmin)sedheadmin=sedhead[i]; 6628 } 6629 6630 for(i=0;i<numdof;i++){ 6631 /*If mask was alread one, keep one*/ 6632 if(old_active[i]>0.){ 6633 vec_mask->SetValue(nodes[i]->Sid(),1.,INS_VAL); 6634 } 6635 6636 this->GetHydrologyDCInefficientHmax(&h_max,nodes[i]); 6637 /*Increase the size of the efficient system if needed*/ 6638 /*Increase is needed if the epl head reach the maximum value (sediment value for now)*/ 6639 if(eplhead[i]>=h_max){ 6640 vec_mask->SetValue(nodes[i]->Sid(),1.,INS_VAL); 6641 for(j=0;j<numdof;j++){ 6642 6643 /*Increase of the domain is on the downstream node in term of sediment head*/ 6644 if(sedhead[j] == sedheadmin){ 6645 vec_mask->SetValue(nodes[j]->Sid(),1.,INS_VAL); 6646 break; 6647 } 6648 } 6652 6649 } 6653 6650 } -
issm/trunk-jpl/src/c/classes/Elements/Tria.h
r15230 r15249 89 89 int Sid(); 90 90 bool IsOnBed(); 91 bool IsFloating(); 91 bool IsFloating(); 92 92 bool IsNodeOnShelf(); 93 93 bool IsNodeOnShelfFromFlags(IssmDouble* flags); … … 196 196 ElementVector* CreatePVectorSlope(void); 197 197 IssmDouble GetArea(void); 198 void BasisIntegral(Vector<IssmDouble>* gbasis);199 198 int GetElementType(void); 200 199 void GetDofList(int** pdoflist,int approximation_enum,int setenum); … … 250 249 void GetSolutionFromInputsHydrologyDCEfficient(Vector<IssmDouble>* solution); 251 250 void CreateHydrologyWaterVelocityInput(void); 252 void HydrologyUpdateEplDomain(void);253 251 void UpdateConstraints(void); 254 252 void InputUpdateFromSolutionHydrology(IssmDouble* solution); … … 259 257 void GetHydrologyDCInefficientHmax(IssmDouble* ph_max, Node* innode); 260 258 void GetHydrologyTransfer(Vector<IssmDouble>* transfer); 259 void HydrologyEPLGetActive(Vector<IssmDouble>* active_vec); 260 void HydrologyEPLGetMask(Vector<IssmDouble>* vec_mask); 261 261 bool AllActive(void); 262 bool AnyActive(void); 262 263 #endif 263 264 #ifdef _HAVE_BALANCED_ -
issm/trunk-jpl/src/c/classes/FemModel.cpp
r15237 r15249 1581 1581 #endif 1582 1582 1583 void FemModel::BasisIntegralsx(void){ /*{{{*/ 1584 1585 Vector<IssmDouble>* basisg=NULL; 1586 1587 /*Vector allocation*/ 1588 basisg=new Vector<IssmDouble>(nodes->NumberOfNodes()); 1589 1583 void FemModel::HydrologyEPLupdateDomainx(void){ /*{{{*/ 1584 1585 Vector<IssmDouble> *mask = NULL; 1586 IssmDouble *serial_mask = NULL; 1587 Vector<IssmDouble> *active = NULL; 1588 IssmDouble *serial_active = NULL; 1589 1590 /*Step 1: update maks, the mask might be extended by residual and/or using downstream sediment head*/ 1591 mask=new Vector<IssmDouble>(nodes->NumberOfNodes()); 1590 1592 for (int i=0;i<elements->Size();i++){ 1591 1593 Element* element=dynamic_cast<Element*>(elements->GetObjectByOffset(i)); 1592 element->BasisIntegral(basisg); 1593 } 1594 /*Assemble*/ 1595 basisg->Assemble(); 1596 1597 /*Update Inputs*/ 1598 InputUpdateFromVectorx(elements,nodes,vertices,loads,materials,parameters,basisg,BasisIntegralEnum,NodesEnum); 1594 element->HydrologyEPLGetMask(mask); 1595 } 1596 1597 /*Assemble and serialize*/ 1598 mask->Assemble(); 1599 serial_mask=mask->ToMPISerial(); 1600 delete mask; 1601 1602 /*Update Mask*/ 1603 InputUpdateFromVectorx(elements,nodes,vertices,loads,materials,parameters,serial_mask,HydrologydcMaskEplactiveEnum,NodeSIdEnum); 1604 xDelete<IssmDouble>(serial_mask); 1605 1606 /*Step 2: update node activity. If one element is connected to mask=1, all nodes are active*/ 1607 active=new Vector<IssmDouble>(nodes->NumberOfNodes()); 1608 for (int i=0;i<elements->Size();i++){ 1609 Element* element=dynamic_cast<Element*>(elements->GetObjectByOffset(i)); 1610 element->HydrologyEPLGetActive(active); 1611 } 1612 1613 /*Assemble and serialize*/ 1614 active->Assemble(); 1615 serial_active=active->ToMPISerial(); 1616 delete active; 1617 1618 /*Update node activation accordingly*/ 1619 for (int i=0;i<nodes->Size();i++){ 1620 Node* node=dynamic_cast<Node*>(nodes->GetObjectByOffset(i)); 1621 if(serial_active[node->Sid()]==1.){ 1622 node->Activate(); 1623 } 1624 else{ 1625 node->Deactivate(); 1626 } 1627 } 1628 xDelete<IssmDouble>(serial_active); 1599 1629 1600 1630 } -
issm/trunk-jpl/src/c/classes/FemModel.h
r15130 r15249 99 99 void UpdateConstraintsx(void); 100 100 int UpdateVertexPositionsx(void); 101 void BasisIntegralsx(void);101 void ParEplMask(void); 102 102 void HydrologyTransferx(void); 103 void HydrologyEPLupdateDomainx(void); 103 104 }; 104 105 -
issm/trunk-jpl/src/c/solutionsequences/solutionsequence_hydro_nonlinear.cpp
r15234 r15249 41 41 femmodel->parameters->FindParam(&eps_hyd,HydrologydcRelTolEnum); 42 42 femmodel->parameters->FindParam(&time,TimeEnum); 43 femmodel->BasisIntegralsx();44 43 hydro_maxiter=100; 45 44 hydrocount=1; … … 123 122 delete uf; 124 123 Solverx(&uf, Kff, pf,old_uf, df, femmodel->parameters); 125 delete old_uf; old_uf=uf->Duplicate(); 124 delete old_uf; old_uf=uf->Duplicate(); 126 125 if(eplcount>1) delete ug_epl; 127 126 delete Kff;delete pf; … … 129 128 Mergesolutionfromftogx(&ug_epl,uf,ys,femmodel->nodes,femmodel->parameters); delete ys; 130 129 InputUpdateFromSolutionx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,ug_epl); 130 femmodel->HydrologyEPLupdateDomainx(); 131 131 ConstraintsStatex(&constraints_converged,&num_unstable_constraints, femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters); 132 132 … … 144 144 InputUpdateFromConstantx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,sediment_kmax,MeltingOffsetEnum); 145 145 InputUpdateFromSolutionx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,ug_epl); 146 femmodel->HydrologyEPLupdateDomainx(); 146 147 break; 147 148 }
Note:
See TracChangeset
for help on using the changeset viewer.