Index: /issm/trunk-jpl/src/c/classes/objects/Elements/Penta.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/objects/Elements/Penta.cpp (revision 14590) +++ /issm/trunk-jpl/src/c/classes/objects/Elements/Penta.cpp (revision 14591) @@ -3487,5 +3487,4 @@ GetNodalFunctionsP1(&L[0], gauss); D_scalar_trans=gauss->weight*Jdet; - D_scalar_trans=D_scalar_trans; TripleMultiply(&L[0],numdof,1,0, @@ -3717,5 +3716,4 @@ GetNodalFunctionsP1(&L[0], gauss); D_scalar_trans=gauss->weight*Jdet; - D_scalar_trans=D_scalar_trans; TripleMultiply(&L[0],numdof,1,0, Index: /issm/trunk-jpl/src/c/classes/objects/Elements/Tria.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/objects/Elements/Tria.cpp (revision 14590) +++ /issm/trunk-jpl/src/c/classes/objects/Elements/Tria.cpp (revision 14591) @@ -533,5 +533,5 @@ IssmDouble D,Jdet; IssmDouble xyz_list[NUMVERTICES][3]; - IssmDouble L[1][3]; + IssmDouble L[1][numdof]; GaussTria *gauss = NULL; @@ -5677,4 +5677,21 @@ ElementMatrix* Tria::CreateKMatrixHydrology(void){ + int hydrology_model; + parameters->FindParam(&hydrology_model,HydrologyEnum); + + switch(hydrology_model){ + case HydrologyshreveEnum: + return CreateKMatrixHydrologyShreve(); + case HydrologydcEnum: + return CreateKMatrixHydrologyDC(); + default: + _error_("Approximation " << EnumToStringx(hydrology_model) << " not supported yet"); + } + +} +/*}}}*/ +/*FUNCTION Tria::CreateKMatrixHydrologyShreve{{{*/ +ElementMatrix* Tria::CreateKMatrixHydrologyShreve(void){ + /*Constants*/ const int numdof=NDOF1*NUMVERTICES; @@ -5778,6 +5795,82 @@ } /*}}}*/ -/*FUNCTION Tria::CreatePVectorHydrology {{{*/ +/*FUNCTION Tria::CreateKMatrixHydrologyDC{{{*/ +ElementMatrix* Tria::CreateKMatrixHydrologyDC(void){ + + /*constants: */ + const int numdof=NDOF1*NUMVERTICES; + + /* Intermediaries */ + IssmDouble D_scalar,Jdet; + IssmDouble sediment_porosity,dt; + IssmDouble xyz_list[NUMVERTICES][3]; + IssmDouble B[2][numdof]; + IssmDouble L[NUMVERTICES]; + IssmDouble D[2][2]; + GaussTria *gauss = NULL; + + /*Initialize Element matrix*/ + ElementMatrix* Ke=new ElementMatrix(nodes,NUMVERTICES,this->parameters,NoneApproximationEnum); + + /*Retrieve all inputs and parameters*/ + GetVerticesCoordinates(&xyz_list[0][0],nodes,NUMVERTICES); + sediment_porosity = matpar->GetSedimentPorosity(); + this->parameters->FindParam(&dt,TimesteppingTimeStepEnum); + + /* 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); + + /*Diffusivity*/ + D_scalar=sediment_porosity*gauss->weight*Jdet; + if(reCast(dt)) D_scalar=D_scalar*dt; + D[0][0]=D_scalar; D[0][1]=0.; + D[1][0]=0.; D[1][1]=D_scalar; + GetBHydro(&B[0][0],&xyz_list[0][0],gauss); + TripleMultiply(&B[0][0],2,numdof,1, + &D[0][0],2,2,0, + &B[0][0],2,numdof,0, + &Ke->values[0],1); + + /*Transient*/ + if(reCast(dt)){ + GetL(&L[0], &xyz_list[0][0], gauss,NDOF1); + D_scalar=gauss->weight*Jdet; + + TripleMultiply(&L[0],numdof,1,0, + &D_scalar,1,1,0, + &L[0],1,numdof,0, + &Ke->values[0],1); + } + } + + /*Clean up and return*/ + delete gauss; + return Ke; +} +/*}}}*/ +/*FUNCTION Tria::CreatePVectorHydrology{{{*/ ElementVector* Tria::CreatePVectorHydrology(void){ + + int hydrology_model; + parameters->FindParam(&hydrology_model,HydrologyEnum); + + switch(hydrology_model){ + case HydrologyshreveEnum: + return CreatePVectorHydrologyShreve(); + case HydrologydcEnum: + return CreatePVectorHydrologyDC(); + default: + _error_("Approximation " << EnumToStringx(hydrology_model) << " not supported yet"); + } + +} +/*}}}*/ +/*FUNCTION Tria::CreatePVectorHydrologyShreve {{{*/ +ElementVector* Tria::CreatePVectorHydrologyShreve(void){ /*Constants*/ @@ -5820,4 +5913,55 @@ if(reCast(dt))for(i=0;ivalues[i]+=Jdettria*gauss->weight*(old_watercolumn_g+dt*basal_melting_g)*basis[i]; else for(i=0;ivalues[i]+=Jdettria*gauss->weight*basal_melting_g*basis[i]; + } + + /*Clean up and return*/ + delete gauss; + return pe; +} +/*}}}*/ +/*FUNCTION Tria::CreatePVectorHydrologyDC {{{*/ +ElementVector* Tria::CreatePVectorHydrologyDC(void){ + + /*Constants*/ + const int numdof=NDOF1*NUMVERTICES; + + /*Intermediaries */ + IssmDouble Jdet; + IssmDouble xyz_list[NUMVERTICES][3]; + IssmDouble dt,gravity,scalar,water_head; + IssmDouble basis[numdof]; + GaussTria* gauss=NULL; + + /*Initialize Element vector*/ + ElementVector* pe=new ElementVector(nodes,NUMVERTICES,this->parameters); + + /*Retrieve all inputs and parameters*/ + GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES); + gravity = matpar->GetG(); + Input* old_wh_input=NULL; + if(reCast(dt)){ + old_wh_input=inputs->GetInput(SedimentHeadEnum); _assert_(old_wh_input); + } + + /* 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, gauss); + + /*Gravity term*/ + scalar = Jdet*gauss->weight*gravity; + if(reCast(dt)) scalar = scalar*dt; + for(int i=0;ivalues[i]+=scalar*basis[i]; + + /*Transient term*/ + if(reCast(dt)){ + old_wh_input->GetInputValue(&water_head,gauss); + scalar = Jdet*gauss->weight*water_head; + for(int i=0;ivalues[i]+=scalar*basis[i]; + } } Index: /issm/trunk-jpl/src/c/classes/objects/Elements/Tria.h =================================================================== --- /issm/trunk-jpl/src/c/classes/objects/Elements/Tria.h (revision 14590) +++ /issm/trunk-jpl/src/c/classes/objects/Elements/Tria.h (revision 14591) @@ -239,5 +239,9 @@ #ifdef _HAVE_HYDROLOGY_ ElementMatrix* CreateKMatrixHydrology(void); + ElementMatrix* CreateKMatrixHydrologyShreve(void); + ElementMatrix* CreateKMatrixHydrologyDC(void); ElementVector* CreatePVectorHydrology(void); + ElementVector* CreatePVectorHydrologyShreve(void); + ElementVector* CreatePVectorHydrologyDC(void); void CreateHydrologyWaterVelocityInput(void); void GetSolutionFromInputsHydrology(Vector* solution); Index: /issm/trunk-jpl/src/c/classes/objects/Elements/TriaRef.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/objects/Elements/TriaRef.cpp (revision 14590) +++ /issm/trunk-jpl/src/c/classes/objects/Elements/TriaRef.cpp (revision 14591) @@ -56,4 +56,29 @@ /*Reference Element numerics*/ +/*FUNCTION TriaRef::GetBHydro {{{*/ +void TriaRef::GetBHydro(IssmDouble* B, IssmDouble* xyz_list, GaussTria* gauss){ + /*Compute B matrix. B=[B1 B2 B3] where Bi is of size 3*NDOF2. + * For node i, Bi can be expressed in the actual coordinate system + * by: + * Bi=[ dh/dx ] + * [ dh/dy ] + * where h is the interpolation function for node i. + * + * We assume B has been allocated already, of size: 3x(NDOF2*NUMNODES) + */ + + int i; + IssmDouble dbasis[NDOF2][NUMNODES]; + + /*Get dh1dh2dh3 in actual coordinate system: */ + GetNodalFunctionsDerivatives(&dbasis[0][0],xyz_list,gauss); + + /*Build B: */ + for (i=0;iSetCurrentConfiguration(HydrologyAnalysisEnum); @@ -80,6 +79,7 @@ } else if (hydrology_model==HydrologydcEnum){ - //solver_linear(femmodel,modify_loads); - _error_("not implemented yet"); + if(VerboseSolution()) _pprintLine_(" computing water head"); + femmodel->SetCurrentConfiguration(HydrologyAnalysisEnum); + solver_linear(femmodel); if(save_results && ((i+1)%output_frequency==0 || (i+1)==nsteps)){ if(VerboseSolution()) _pprintLine_(" saving results ");