Index: /issm/trunk-jpl/src/c/analyses/HydrologyShreveAnalysis.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/HydrologyShreveAnalysis.cpp (revision 16902) +++ /issm/trunk-jpl/src/c/analyses/HydrologyShreveAnalysis.cpp (revision 16903) @@ -83,5 +83,92 @@ /*Finite Element Analysis*/ ElementMatrix* HydrologyShreveAnalysis::CreateKMatrix(Element* element){/*{{{*/ - _error_("not implemented yet"); + + /*Intermediaries */ + IssmDouble diffusivity; + IssmDouble Jdet,D_scalar,dt,h; + IssmDouble vx,vy,vel,dvxdx,dvydy; + IssmDouble dvx[2],dvy[2]; + IssmDouble* xyz_list = NULL; + + /*Fetch number of nodes and dof for this finite element*/ + int numnodes = element->GetNumberOfNodes(); + + /*Initialize Element vector and other vectors*/ + ElementMatrix* Ke = element->NewElementMatrix(); + IssmDouble* basis = xNew(numnodes); + IssmDouble* B = xNew(2*numnodes); + IssmDouble* Bprime = xNew(2*numnodes); + IssmDouble D[2][2]={0.}; + + /*Create water velocity vx and vy from current inputs*/ + CreateHydrologyWaterVelocityInput(element); + + /*Retrieve all inputs and parameters*/ + element->GetVerticesCoordinates(&xyz_list); + element->FindParam(&dt,TimesteppingTimeStepEnum); + Input* vx_input=element->GetInput(HydrologyWaterVxEnum); _assert_(vx_input); + Input* vy_input=element->GetInput(HydrologyWaterVyEnum); _assert_(vy_input); + h = element->CharacteristicLength(); + + /* Start looping on the number of gaussian points: */ + Gauss* gauss=element->NewGauss(2); + for(int ig=gauss->begin();igend();ig++){ + gauss->GaussPoint(ig); + + element->JacobianDeterminant(&Jdet,xyz_list,gauss); + element->NodalFunctions(basis,gauss); + + vx_input->GetInputValue(&vx,gauss); + vy_input->GetInputValue(&vy,gauss); + vx_input->GetInputDerivativeValue(&dvx[0],xyz_list,gauss); + vy_input->GetInputDerivativeValue(&dvy[0],xyz_list,gauss); + + D_scalar=gauss->weight*Jdet; + + TripleMultiply(basis,1,numnodes,1, + &D_scalar,1,1,0, + basis,1,numnodes,0, + Ke->values,1); + + GetB(B,element,xyz_list,gauss); + GetBprime(Bprime,element,xyz_list,gauss); + + dvxdx=dvx[0]; + dvydy=dvy[1]; + D_scalar=dt*gauss->weight*Jdet; + + D[0][0]=D_scalar*dvxdx; + D[1][1]=D_scalar*dvydy; + TripleMultiply(B,2,numnodes,1, + &D[0][0],2,2,0, + B,2,numnodes,0, + &Ke->values[0],1); + + D[0][0]=D_scalar*vx; + D[1][1]=D_scalar*vy; + TripleMultiply(B,2,numnodes,1, + &D[0][0],2,2,0, + Bprime,2,numnodes,0, + &Ke->values[0],1); + + /*Artificial diffusivity*/ + vel=sqrt(vx*vx+vy*vy); + D[0][0]=D_scalar*diffusivity*h/(2*vel)*vx*vx; + D[1][0]=D_scalar*diffusivity*h/(2*vel)*vy*vx; + D[0][1]=D_scalar*diffusivity*h/(2*vel)*vx*vy; + D[1][1]=D_scalar*diffusivity*h/(2*vel)*vy*vy; + TripleMultiply(Bprime,2,numnodes,1, + &D[0][0],2,2,0, + Bprime,2,numnodes,0, + &Ke->values[0],1); + } + + /*Clean up and return*/ + xDelete(xyz_list); + xDelete(basis); + xDelete(B); + xDelete(Bprime); + delete gauss; + return Ke; }/*}}}*/ ElementVector* HydrologyShreveAnalysis::CreatePVector(Element* element){/*{{{*/ @@ -134,4 +221,59 @@ return pe; }/*}}}*/ +void HydrologyShreveAnalysis::GetB(IssmDouble* B,Element* element,IssmDouble* xyz_list,Gauss* 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=[ N ] + * [ N ] + * where N is the finiteelement function for node i. + * + * We assume B_prog has been allocated already, of size: 2x(NDOF1*numnodes) + */ + + /*Fetch number of nodes for this finite element*/ + int numnodes = element->GetNumberOfNodes(); + + /*Get nodal functions*/ + IssmDouble* basis=xNew(numnodes); + element->NodalFunctions(basis,gauss); + + /*Build B: */ + for(int i=0;i(basis); +}/*}}}*/ +void HydrologyShreveAnalysis::GetBprime(IssmDouble* Bprime,Element* element,IssmDouble* xyz_list,Gauss* 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_prime=[ dN/dx ] + * [ dN/dy ] + * where N is the finiteelement function for node i. + * + * We assume B' has been allocated already, of size: 3x(NDOF2*numnodes) + */ + + /*Fetch number of nodes for this finite element*/ + int numnodes = element->GetNumberOfNodes(); + + /*Get nodal functions derivatives*/ + IssmDouble* dbasis=xNew(2*numnodes); + element->NodalFunctionsDerivatives(dbasis,xyz_list,gauss); + + /*Build B': */ + for(int i=0;i(dbasis); + +}/*}}}*/ void HydrologyShreveAnalysis::GetSolutionFromInputs(Vector* solution,Element* element){/*{{{*/ element->GetSolutionFromInputsOneDof(solution,WatercolumnEnum); @@ -163,2 +305,55 @@ xDelete(doflist); }/*}}}*/ + +/*Intermediaries*/ +void HydrologyShreveAnalysis::CreateHydrologyWaterVelocityInput(Element* element){/*{{{*/ + + /*Intermediaries*/ + IssmDouble dsdx,dsdy,dbdx,dbdy,w; + + /*Retrieve all inputs and parameters*/ + IssmDouble rho_ice = element->GetMaterialParameter(MaterialsRhoIceEnum); + IssmDouble rho_water = element->GetMaterialParameter(MaterialsRhoWaterEnum); + IssmDouble g = element->GetMaterialParameter(ConstantsGEnum); + IssmDouble CR = element->GetMaterialParameter(HydrologyshreveCREnum); + IssmDouble n_man = element->GetMaterialParameter(HydrologyshreveNEnum); + IssmDouble mu_water = element->GetMaterialParameter(MaterialsMuWaterEnum); + Input* surfaceslopex_input = element->GetInput(SurfaceSlopeXEnum); _assert_(surfaceslopex_input); + Input* surfaceslopey_input = element->GetInput(SurfaceSlopeYEnum); _assert_(surfaceslopey_input); + Input* bedslopex_input = element->GetInput(BedSlopeXEnum); _assert_(bedslopex_input); + Input* bedslopey_input = element->GetInput(BedSlopeYEnum); _assert_(bedslopey_input); + Input* watercolumn_input = element->GetInput(WatercolumnEnum); _assert_(watercolumn_input); + + /* compute VelocityFactor */ + IssmDouble VelocityFactor = n_man*CR*CR*rho_water*g/mu_water; + + /*Fetch number of vertices and allocate output*/ + int numvertices = element->GetNumberOfVertices(); + IssmDouble* vx = xNew(numvertices); + IssmDouble* vy = xNew(numvertices); + + Gauss* gauss=element->NewGauss(); + for(int iv=0;ivGaussVertex(iv); + surfaceslopex_input->GetInputValue(&dsdx,gauss); + surfaceslopey_input->GetInputValue(&dsdy,gauss); + bedslopex_input->GetInputValue(&dbdx,gauss); + bedslopey_input->GetInputValue(&dbdy,gauss); + watercolumn_input->GetInputValue(&w,gauss); + + /* Water velocity x and y components */ + // vx[iv]= - w*w/(12 * mu_water)*(rho_ice*g*dsdx+(rho_water-rho_ice)*g*dbdx); + // vy[iv]= - w*w/(12 * mu_water)*(rho_ice*g*dsdy+(rho_water-rho_ice)*g*dbdy); + vx[iv]= - w*w/(VelocityFactor* mu_water)*(rho_ice*g*dsdx+(rho_water-rho_ice)*g*dbdx); + vy[iv]= - w*w/(VelocityFactor* mu_water)*(rho_ice*g*dsdy+(rho_water-rho_ice)*g*dbdy); + } + + /*clean-up*/ + delete gauss; + + /*Add to inputs*/ + element->AddInput(HydrologyWaterVxEnum,vx,P1Enum); + element->AddInput(HydrologyWaterVyEnum,vy,P1Enum); + xDelete(vx); + xDelete(vy); +}/*}}}*/ Index: /issm/trunk-jpl/src/c/analyses/HydrologyShreveAnalysis.h =================================================================== --- /issm/trunk-jpl/src/c/analyses/HydrologyShreveAnalysis.h (revision 16902) +++ /issm/trunk-jpl/src/c/analyses/HydrologyShreveAnalysis.h (revision 16903) @@ -25,4 +25,9 @@ void GetSolutionFromInputs(Vector* solution,Element* element); void InputUpdateFromSolution(IssmDouble* solution,Element* element); + + /*Intermediaries*/ + void GetB(IssmDouble* B,Element* element,IssmDouble* xyz_list,Gauss* gauss); + void GetBprime(IssmDouble* Bprime,Element* element,IssmDouble* xyz_list,Gauss* gauss); + void CreateHydrologyWaterVelocityInput(Element* element); }; #endif Index: /issm/trunk-jpl/src/c/analyses/MasstransportAnalysis.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/MasstransportAnalysis.cpp (revision 16902) +++ /issm/trunk-jpl/src/c/analyses/MasstransportAnalysis.cpp (revision 16903) @@ -252,5 +252,5 @@ IssmDouble* B = xNew(2*numnodes); IssmDouble* Bprime = xNew(2*numnodes); - IssmDouble D[2][2]; + IssmDouble D[2][2]={0.}; /*Retrieve all inputs and parameters*/ @@ -298,6 +298,4 @@ D[0][0]=D_scalar*dvxdx; - D[0][1]=0.; - D[1][0]=0.; D[1][1]=D_scalar*dvydy; TripleMultiply(B,2,numnodes,1, @@ -326,6 +324,4 @@ vyaverage_input->GetInputAverage(&vy); D[0][0]=h/2.0*fabs(vx); - D[0][1]=0.; - D[1][0]=0.; D[1][1]=h/2.0*fabs(vy); } Index: /issm/trunk-jpl/src/c/classes/Materials/Matpar.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Materials/Matpar.cpp (revision 16902) +++ /issm/trunk-jpl/src/c/classes/Materials/Matpar.cpp (revision 16903) @@ -241,25 +241,30 @@ switch(enum_in){ - case MaterialsRhoIceEnum: return this->rho_ice; - case MaterialsRhoWaterEnum: return this->rho_water; - case MaterialsRhoFreshwaterEnum: return this->rho_freshwater; - case MaterialsMuWaterEnum: return this->mu_water; - case MaterialsHeatcapacityEnum: return this->heatcapacity; - case MaterialsThermalconductivityEnum: return this->thermalconductivity; - case MaterialsTemperateiceconductivityEnum: return this->temperateiceconductivity; - case MaterialsLatentheatEnum: return this->latentheat; - case MaterialsBetaEnum: return this->beta; - case MaterialsMeltingpointEnum: return this->meltingpoint; - case ConstantsReferencetemperatureEnum: return this->referencetemperature; - case MaterialsMixedLayerCapacityEnum: return this->mixed_layer_capacity; - case MaterialsThermalExchangeVelocityEnum: return this->thermal_exchange_velocity; - case HydrologydcSedimentPorosityEnum: return this->sediment_porosity; - case HydrologydcSedimentThicknessEnum: return this->sediment_thickness; - case HydrologydcSedimentCompressibilityEnum:return this->sediment_compressibility; - case HydrologydcEplPorosityEnum: return this->epl_porosity; - case HydrologydcEplCompressibilityEnum: return this->epl_compressibility; - case HydrologydcWaterCompressibilityEnum: return this->water_compressibility; - case HydrologydcSedimentTransmitivityEnum: return this->sediment_transmitivity; - case ConstantsGEnum: return this->g; + case MaterialsRhoIceEnum: return this->rho_ice; + case MaterialsRhoWaterEnum: return this->rho_water; + case MaterialsRhoFreshwaterEnum: return this->rho_freshwater; + case MaterialsMuWaterEnum: return this->mu_water; + case MaterialsHeatcapacityEnum: return this->heatcapacity; + case MaterialsThermalconductivityEnum: return this->thermalconductivity; + case MaterialsTemperateiceconductivityEnum: return this->temperateiceconductivity; + case MaterialsLatentheatEnum: return this->latentheat; + case MaterialsBetaEnum: return this->beta; + case MaterialsMeltingpointEnum: return this->meltingpoint; + case ConstantsReferencetemperatureEnum: return this->referencetemperature; + case MaterialsMixedLayerCapacityEnum: return this->mixed_layer_capacity; + case MaterialsThermalExchangeVelocityEnum: return this->thermal_exchange_velocity; + case HydrologydcSedimentPorosityEnum: return this->sediment_porosity; + case HydrologydcSedimentThicknessEnum: return this->sediment_thickness; + case HydrologydcSedimentCompressibilityEnum: return this->sediment_compressibility; + case HydrologydcEplPorosityEnum: return this->epl_porosity; + case HydrologydcEplCompressibilityEnum: return this->epl_compressibility; + case HydrologydcWaterCompressibilityEnum: return this->water_compressibility; + case HydrologydcSedimentTransmitivityEnum: return this->sediment_transmitivity; + case HydrologyshreveCREnum: return this->hydro_CR; + case HydrologyshreveKnEnum: return this->hydro_kn; + case HydrologyshreveNEnum: return this->hydro_n; + case HydrologyshrevePEnum: return this->hydro_p; + case HydrologyshreveQEnum: return this->hydro_q; + case ConstantsGEnum: return this->g; default: _error_("Enum "<