Index: /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp (revision 17241) +++ /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp (revision 17242) @@ -18,5 +18,11 @@ case SSAApproximationEnum: numdofs =2; break; case L1L2ApproximationEnum: numdofs =2; break; - case HOApproximationEnum: numdofs =2; break; + case HOApproximationEnum: + switch(meshtype){ + case Mesh3DEnum: numdofs=2; break; + case Mesh2DverticalEnum: numdofs=1; break; + default: _error_("mesh type not supported yet"); + } + break; case SIAApproximationEnum: numdofs =2; break; case FSvelocityEnum: @@ -478,5 +484,7 @@ else{ IoModelToConstraintsx(constraints,iomodel,StressbalanceSpcvxEnum,StressbalanceAnalysisEnum,finiteelement,1); - IoModelToConstraintsx(constraints,iomodel,StressbalanceSpcvyEnum,StressbalanceAnalysisEnum,finiteelement,2); + if(iomodel->meshtype==Mesh3DEnum){ + IoModelToConstraintsx(constraints,iomodel,StressbalanceSpcvyEnum,StressbalanceAnalysisEnum,finiteelement,2); + } } @@ -955,10 +963,18 @@ IssmDouble vx,vy; - int approximation; + int meshtype,dim,approximation; int* doflist = NULL; + + /*Get some parameters*/ + element->FindParam(&meshtype,MeshTypeEnum); + switch(meshtype){ + case Mesh2DverticalEnum: dim = 2; break; + case Mesh3DEnum: dim = 3; break; + default: _error_("mesh "<GetNumberOfNodes(); - int numdof = numnodes*2; + int numdof = numnodes*(dim-1); element->GetInputValue(&approximation,ApproximationEnum); @@ -969,5 +985,6 @@ /*Get inputs*/ Input* vx_input=element->GetInput(VxEnum); _assert_(vx_input); - Input* vy_input=element->GetInput(VyEnum); _assert_(vy_input); + Input* vy_input=NULL; + if(dim==3){vy_input=element->GetInput(VyEnum); _assert_(vy_input);} /*Ok, we have vx and vy in values, fill in vx and vy arrays: */ @@ -978,7 +995,9 @@ /*Recover vx and vy*/ vx_input->GetInputValue(&vx,gauss); - vy_input->GetInputValue(&vy,gauss); - values[i*2+0]=vx; - values[i*2+1]=vy; + values[i*(dim-1)+0]=vx; + if(dim==3){ + vy_input->GetInputValue(&vy,gauss); + values[i*(dim-1)+1]=vy; + } } @@ -2021,4 +2040,5 @@ /*Intermediaries*/ + int dim,meshtype,bsize; IssmDouble viscosity,newviscosity,oldviscosity; IssmDouble viscosity_overshoot,thickness,Jdet; @@ -2026,4 +2046,12 @@ IssmDouble *xyz_list = NULL; + /*Get problem dimension*/ + element->FindParam(&meshtype,MeshTypeEnum); + switch(meshtype){ + case Mesh2DverticalEnum: dim = 2; bsize = 2; break; + case Mesh3DEnum: dim = 3; bsize = 5; break; + default: _error_("mesh "<GetNumberOfNodes(); @@ -2032,14 +2060,18 @@ /*Initialize Element matrix and vectors*/ ElementMatrix* Ke = element->NewElementMatrix(HOApproximationEnum); - IssmDouble* B = xNew(5*numdof); - IssmDouble* Bprime = xNew(5*numdof); - IssmDouble* D = xNewZeroInit(5*5); + IssmDouble* B = xNew(bsize*numdof); + IssmDouble* Bprime = xNew(bsize*numdof); + IssmDouble* D = xNewZeroInit(bsize*bsize); /*Retrieve all inputs and parameters*/ element->GetVerticesCoordinates(&xyz_list); - Input* vx_input=element->GetInput(VxEnum); _assert_(vx_input); - Input* vy_input=element->GetInput(VyEnum); _assert_(vy_input); - Input* vxold_input=element->GetInput(VxPicardEnum); _assert_(vxold_input); - Input* vyold_input=element->GetInput(VyPicardEnum); _assert_(vyold_input); + Input* vx_input = element->GetInput(VxEnum); _assert_(vx_input); + Input* vxold_input = element->GetInput(VxPicardEnum); _assert_(vxold_input); + Input* vy_input = NULL; + Input* vyold_input = NULL; + if(dim==3){ + vy_input=element->GetInput(VyEnum); _assert_(vy_input); + vyold_input=element->GetInput(VyPicardEnum); _assert_(vyold_input); + } element->FindParam(&viscosity_overshoot,StressbalanceViscosityOvershootEnum); @@ -2050,17 +2082,17 @@ element->JacobianDeterminant(&Jdet,xyz_list,gauss); - this->GetBHO(B,element,xyz_list,gauss); - this->GetBHOprime(Bprime,element,xyz_list,gauss); - - element->ViscosityHO(&viscosity,xyz_list,gauss,vx_input,vy_input); - element->ViscosityHO(&oldviscosity,xyz_list,gauss,vxold_input,vyold_input); + this->GetBHO(B,element,dim,xyz_list,gauss); + this->GetBHOprime(Bprime,element,dim,xyz_list,gauss); + + element->ViscosityHO(&viscosity,dim,xyz_list,gauss,vx_input,vy_input); + element->ViscosityHO(&oldviscosity,dim,xyz_list,gauss,vxold_input,vyold_input); newviscosity=viscosity+viscosity_overshoot*(viscosity-oldviscosity); D_scalar=2.*newviscosity*gauss->weight*Jdet; - for(int i=0;i<5;i++) D[i*5+i]=D_scalar; - - TripleMultiply(B,5,numdof,1, - D,5,5,0, - Bprime,5,numdof,0, + for(int i=0;ivalues[0],1); } @@ -2082,4 +2114,5 @@ /*Intermediaries*/ + int dim,meshtype; bool mainlyfloating; int migration_style,point1; @@ -2089,20 +2122,29 @@ Gauss* gauss = NULL; + /*Get problem dimension*/ + element->FindParam(&meshtype,MeshTypeEnum); + switch(meshtype){ + case Mesh2DverticalEnum: dim = 2; break; + case Mesh3DEnum: dim = 3; break; + default: _error_("mesh "<GetNumberOfNodes(); - int numdof = numnodes*2; + int numdof = numnodes*(dim-1); /*Initialize Element matrix and vectors*/ - ElementMatrix* Ke = element->NewElementMatrix(HOApproximationEnum); - IssmDouble* B = xNew(2*numdof); - IssmDouble D[2][2] = {0.}; + ElementMatrix* Ke = element->NewElementMatrix(HOApproximationEnum); + IssmDouble* B = xNew((dim-1)*numdof); + IssmDouble* D = xNew((dim-1)*(dim-1)); /*Retrieve all inputs and parameters*/ element->GetVerticesCoordinatesBase(&xyz_list_base); element->FindParam(&migration_style,GroundinglineMigrationEnum); - Input* vx_input = element->GetInput(VxEnum); _assert_(vx_input); - Input* vy_input = element->GetInput(VyEnum); _assert_(vy_input); - Input* vz_input = element->GetInput(VzEnum); _assert_(vz_input); + Input* vx_input = element->GetInput(VxEnum); _assert_(vx_input); + Input* vy_input = element->GetInput(VyEnum); _assert_(vy_input); + Input* vz_input = NULL; Input* gllevelset_input = NULL; + if(dim==3){vz_input=element->GetInput(VzEnum); _assert_(vz_input);} /*build friction object, used later on: */ @@ -2132,11 +2174,11 @@ } - this->GetBHOFriction(B,element,xyz_list_base,gauss); + this->GetBHOFriction(B,element,dim,xyz_list_base,gauss); element->JacobianDeterminantBase(&Jdet,xyz_list_base,gauss); - for(int i=0;i<2;i++) D[i][i]=alpha2*gauss->weight*Jdet; - - TripleMultiply(B,2,numdof,1, - &D[0][0],2,2,0, - B,2,numdof,0, + for(int i=0;iweight*Jdet; + + TripleMultiply(B,dim-1,numdof,1, + D,dim-1,dim-1,0, + B,dim-1,numdof,0, &Ke->values[0],1); } @@ -2150,4 +2192,5 @@ xDelete(xyz_list_base); xDelete(B); + xDelete(D); return Ke; }/*}}}*/ @@ -2167,6 +2210,15 @@ /*Intermediaries */ + int dim,meshtype; IssmDouble Jdet,slope[3]; IssmDouble* xyz_list = NULL; + + /*Get problem dimension*/ + element->FindParam(&meshtype,MeshTypeEnum); + switch(meshtype){ + case Mesh2DverticalEnum: dim = 2; break; + case Mesh3DEnum: dim = 3; break; + default: _error_("mesh "<values[i*2+0]+=-rhog*slope[0]*Jdet*gauss->weight*basis[i]; - pe->values[i*2+1]+=-rhog*slope[1]*Jdet*gauss->weight*basis[i]; + pe->values[i*(dim-1)+0]+=-rhog*slope[0]*Jdet*gauss->weight*basis[i]; + if(dim==3) pe->values[i*(dim-1)+1]+=-rhog*slope[1]*Jdet*gauss->weight*basis[i]; } } @@ -2212,4 +2264,5 @@ /*Intermediaries*/ + int dim,meshtype; IssmDouble Jdet,surface,z,water_pressure,ice_pressure; IssmDouble surface_under_water,base_under_water,pressure; @@ -2218,4 +2271,12 @@ IssmDouble normal[3]; Gauss* gauss = NULL; + + /*Get problem dimension*/ + element->FindParam(&meshtype,MeshTypeEnum); + switch(meshtype){ + case Mesh2DverticalEnum: dim = 2; break; + case Mesh3DEnum: dim = 3; break; + default: _error_("mesh "<zmax) zmax=xyz_list[i*3+2]; - IssmDouble zmin=xyz_list[0*3+2]; for(int i=1;i0 && zmin<0) gauss=element->NewGauss(xyz_list,xyz_list_front,3,10);//refined in vertical because of the sea level discontinuity - else gauss=element->NewGauss(xyz_list,xyz_list_front,3,3); + IssmDouble zmax=xyz_list[0*3+(dim-1)]; for(int i=1;izmax) zmax=xyz_list[i*3+(dim-1)]; + IssmDouble zmin=xyz_list[0*3+(dim-1)]; for(int i=1;i0. && zmin<0.) gauss=element->NewGauss(xyz_list,xyz_list_front,3,10);//refined in vertical because of the sea level discontinuity + else gauss=element->NewGauss(xyz_list,xyz_list_front,3,3); /* Start looping on the number of gaussian points: */ @@ -2256,6 +2317,6 @@ for (int i=0;ivalues[2*i+0]+= pressure*Jdet*gauss->weight*normal[0]*basis[i]; - pe->values[2*i+1]+= pressure*Jdet*gauss->weight*normal[1]*basis[i]; + pe->values[(dim-1)*i+0]+= pressure*Jdet*gauss->weight*normal[0]*basis[i]; + if(dim==3) pe->values[(dim-1)*i+1]+= pressure*Jdet*gauss->weight*normal[1]*basis[i]; } } @@ -2271,13 +2332,15 @@ return pe; }/*}}}*/ -void StressbalanceAnalysis::GetBHO(IssmDouble* B,Element* element,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/ +void StressbalanceAnalysis::GetBHO(IssmDouble* B,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/ /*Compute B matrix. B=[B1 B2 B3 B4 B5 B6] where Bi is of size 5*NDOF2. * For node i, Bi can be expressed in the actual coordinate system * by: - * Bi=[ dh/dx 0 ] - * [ 0 dh/dy ] - * [ 1/2*dh/dy 1/2*dh/dx ] - * [ 1/2*dh/dz 0 ] - * [ 0 1/2*dh/dz ] + * 3D 2D + * + * Bi=[ dh/dx 0 ] Bi=[ dh/dx] + * [ 0 dh/dy ] [ dh/dy] + * [ 1/2*dh/dy 1/2*dh/dx ] + * [ 1/2*dh/dz 0 ] + * [ 0 1/2*dh/dz ] * where h is the interpolation function for node i. * @@ -2289,31 +2352,41 @@ /*Get nodal functions derivatives*/ - IssmDouble* dbasis=xNew(3*numnodes); + IssmDouble* dbasis=xNew(dim*numnodes); element->NodalFunctionsDerivatives(dbasis,xyz_list,gauss); /*Build B: */ - for(int i=0;i(dbasis); }/*}}}*/ -void StressbalanceAnalysis::GetBHOprime(IssmDouble* Bprime,Element* element,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/ +void StressbalanceAnalysis::GetBHOprime(IssmDouble* Bprime,Element* element,int dim,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=[ 2*dN/dx dN/dy ] - * [ dN/dx 2*dN/dy ] - * [ dN/dy dN/dx ] + * 3D 2D + * Bi_prime=[ 2*dN/dx dN/dy ] Bi_prime=[ 2*dN/dx ] + * [ dN/dx 2*dN/dy ] [ dN/dy ] + * [ dN/dy dN/dx ] * where hNis the finiteelement function for node i. * @@ -2325,29 +2398,38 @@ /*Get nodal functions derivatives*/ - IssmDouble* dbasis=xNew(3*numnodes); + IssmDouble* dbasis=xNew(dim*numnodes); element->NodalFunctionsDerivatives(dbasis,xyz_list,gauss); /*Build B': */ - for(int i=0;i(dbasis); }/*}}}*/ -void StressbalanceAnalysis::GetBHOFriction(IssmDouble* B,Element* element,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/ +void StressbalanceAnalysis::GetBHOFriction(IssmDouble* B,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/ /*Compute B matrix. B=[B1 B2 B3] where Bi is square and of size 2. * For node i, Bi can be expressed in the actual coordinate system * by: - * Bi=[ N 0 ] + * 3D 2D + * Bi=[ N 0 ] Bi=N * [ 0 N ] * where N is the finiteelement function for node i. @@ -2364,9 +2446,16 @@ /*Build L: */ - for(int i=0;iFindParam(&meshtype,MeshTypeEnum); + switch(meshtype){ + case Mesh2DverticalEnum: dim = 2; break; + case Mesh3DEnum: dim = 3; break; + default: _error_("mesh "<GetNumberOfNodes(); - int numdof = numnodes*2; + int numdof = numnodes*(dim-1); /*Fetch dof list and allocate solution vectors*/ element->GetDofList(&doflist,HOApproximationEnum,GsetEnum); - IssmDouble* values = xNew(numdof); - IssmDouble* vx = xNew(numnodes); - IssmDouble* vy = xNew(numnodes); - IssmDouble* vz = xNew(numnodes); - IssmDouble* vel = xNew(numnodes); + IssmDouble* values = xNew(numdof); + IssmDouble* vx = xNew(numnodes); + IssmDouble* vy = xNew(numnodes); + IssmDouble* vz = xNew(numnodes); + IssmDouble* vel = xNew(numnodes); /*Use the dof list to index into the solution vector: */ @@ -2413,21 +2510,26 @@ /*Ok, we have vx and vy in values, fill in vx and vy arrays: */ for(i=0;i(vx[i])) _error_("NaN found in solution vector"); - if(xIsNan(vy[i])) _error_("NaN found in solution vector"); + if(dim==3){ + vy[i]=values[i*(dim-1)+1]; + if(xIsNan(vy[i])) _error_("NaN found in solution vector"); + } } /*Get Vz and compute vel*/ - element->GetInputListOnNodes(&vz[0],VzEnum,0.); - for(i=0;iGetInputListOnNodes(&vz[0],VzEnum,0.); + for(i=0;iGetInputListOnNodes(&vy[0],VyEnum,0.); + for(i=0;iInputChangeName(VxEnum,VxPicardEnum); - element->InputChangeName(VyEnum,VyPicardEnum); - element->InputChangeName(PressureEnum,PressurePicardEnum); + if(dim==3)element->InputChangeName(VyEnum,VyPicardEnum); /*Add vx and vy as inputs to the element: */ @@ -2436,5 +2538,5 @@ //element->AddInput(VelEnum,vel,element->GetElementType()); element->AddInput(VxEnum,vx,P1Enum); - element->AddInput(VyEnum,vy,P1Enum); + if(dim==3)element->AddInput(VyEnum,vy,P1Enum); element->AddInput(VelEnum,vel,P1Enum); @@ -3579,6 +3681,6 @@ if(approximation==SSAHOApproximationEnum){ - element->ViscosityHO(&viscosity,xyz_list,gauss,vx_input,vy_input); - element->ViscosityHO(&oldviscosity,xyz_list,gauss,vxold_input,vyold_input); + element->ViscosityHO(&viscosity,dim,xyz_list,gauss,vx_input,vy_input); + element->ViscosityHO(&oldviscosity,dim,xyz_list,gauss,vxold_input,vyold_input); newviscosity=viscosity+viscosity_overshoot*(viscosity-oldviscosity); } @@ -3793,5 +3895,5 @@ friction->GetAlpha2(&alpha2,gauss,vx_input,vy_input,vz_input); element->JacobianDeterminantBase(&Jdet2d, xyz_list_tria,gauss); - this->GetBHOFriction(L,element,xyz_list_tria,gauss); + this->GetBHOFriction(L,element,3,xyz_list_tria,gauss); DL_scalar=alpha2*gauss->weight*Jdet2d; @@ -3878,6 +3980,6 @@ this->GetBSSAHO(B, element,xyz_list, gauss); this->GetBSSAprime(Bprime, basaltria,xyz_list, gauss_tria); - element->ViscosityHO(&viscosity,xyz_list,gauss,vx_input,vy_input); - element->ViscosityHO(&oldviscosity,xyz_list,gauss,vxold_input,vyold_input); + element->ViscosityHO(&viscosity,3,xyz_list,gauss,vx_input,vy_input); + element->ViscosityHO(&oldviscosity,3,xyz_list,gauss,vxold_input,vyold_input); newviscosity=viscosity+viscosity_overshoot*(viscosity-oldviscosity); Index: /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.h =================================================================== --- /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.h (revision 17241) +++ /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.h (revision 17242) @@ -59,7 +59,7 @@ ElementVector* CreatePVectorHODrivingStress(Element* element); ElementVector* CreatePVectorHOFront(Element* element); - void GetBHO(IssmDouble* B,Element* element,IssmDouble* xyz_list,Gauss* gauss); - void GetBHOprime(IssmDouble* B,Element* element,IssmDouble* xyz_list,Gauss* gauss); - void GetBHOFriction(IssmDouble* B,Element* element,IssmDouble* xyz_list,Gauss* gauss); + void GetBHO(IssmDouble* B,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss); + void GetBHOprime(IssmDouble* B,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss); + void GetBHOFriction(IssmDouble* B,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss); void InputUpdateFromSolutionHO(IssmDouble* solution,Element* element); /*FS*/ Index: /issm/trunk-jpl/src/c/classes/Elements/Element.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Element.cpp (revision 17241) +++ /issm/trunk-jpl/src/c/classes/Elements/Element.cpp (revision 17242) @@ -378,12 +378,19 @@ *pviscosity = viscosity; }/*}}}*/ -void Element::ViscosityHO(IssmDouble* pviscosity,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input){/*{{{*/ +void Element::ViscosityHO(IssmDouble* pviscosity,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input){/*{{{*/ /*Intermediaries*/ IssmDouble viscosity; - IssmDouble epsilon[5];/* epsilon=[exx,eyy,exy,exz,eyz];*/ - - this->StrainRateHO(&epsilon[0],xyz_list,gauss,vx_input,vy_input); - material->GetViscosity3d(&viscosity, &epsilon[0]); + IssmDouble epsilon3d[5];/* epsilon=[exx,eyy,exy,exz,eyz];*/ + IssmDouble epsilon2d[2]; /* epsilon=[exx,exy]; */ + + if(dim==3){ + this->StrainRateHO(&epsilon3d[0],xyz_list,gauss,vx_input,vy_input); + material->GetViscosity3d(&viscosity, &epsilon3d[0]); + } + else{ + this->StrainRateHO2dvertical(&epsilon2d[0],xyz_list,gauss,vx_input,vy_input); + material->GetViscosity2dverticalHO(&viscosity, &epsilon2d[0]); + } /*Assign output pointer*/ @@ -469,4 +476,28 @@ }/*}}}*/ +void Element::StrainRateHO2dvertical(IssmDouble* epsilon,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input){/*{{{*/ + /*Compute the 2d Blatter/HOStrain Rate (2 components): + * + * epsilon=[exx exz] + * + * with exz=1/2 du/dz + * + * the contribution of vz is neglected + */ + + /*Intermediaries*/ + IssmDouble dvx[3]; + + /*Check that both inputs have been found*/ + if (!vx_input){ + _error_("Input missing. Here are the input pointers we have for vx: " << vx_input <<"\n"); + } + + /*Get strain rate assuming that epsilon has been allocated*/ + vx_input->GetInputDerivativeValue(&dvx[0],xyz_list,gauss); + epsilon[0] = dvx[0]; + epsilon[1] = 0.5*dvx[1]; + +}/*}}}*/ void Element::StrainRateSSA(IssmDouble* epsilon,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input){/*{{{*/ Index: /issm/trunk-jpl/src/c/classes/Elements/Element.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Element.h (revision 17241) +++ /issm/trunk-jpl/src/c/classes/Elements/Element.h (revision 17242) @@ -70,9 +70,10 @@ void StrainRateSSA(IssmDouble* epsilon,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input); void StrainRateHO(IssmDouble* epsilon,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input); + void StrainRateHO2dvertical(IssmDouble* epsilon,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input); void StrainRateFS(IssmDouble* epsilon,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input,Input* vz_input); IssmDouble TMeltingPoint(IssmDouble pressure); void ViscousHeatingCreateInput(void); void ViscosityFS(IssmDouble* pviscosity,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input,Input* vz_input); - void ViscosityHO(IssmDouble* pviscosity,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input); + void ViscosityHO(IssmDouble* pviscosity,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input); void ViscosityL1L2(IssmDouble* pviscosity,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input,Input* surf); void ViscositySSA(IssmDouble* pviscosity,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input); Index: /issm/trunk-jpl/src/c/classes/Elements/TriaRef.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/TriaRef.cpp (revision 17241) +++ /issm/trunk-jpl/src/c/classes/Elements/TriaRef.cpp (revision 17242) @@ -451,4 +451,9 @@ switch(finiteelement){ + case P0Enum: + /*Nodal function 1*/ + dbasis[NUMNODESP0*0+0] = 0.; + dbasis[NUMNODESP0*1+0] = 0.; + return; case P1Enum: case P1DGEnum: /*Nodal function 1*/ Index: /issm/trunk-jpl/src/c/classes/Materials/Material.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Materials/Material.h (revision 17241) +++ /issm/trunk-jpl/src/c/classes/Materials/Material.h (revision 17242) @@ -28,4 +28,5 @@ virtual void GetViscosity2d(IssmDouble* pviscosity, IssmDouble* pepsilon)=0; virtual void GetViscosity2dvertical(IssmDouble* pviscosity, IssmDouble* pepsilon)=0; + virtual void GetViscosity2dverticalHO(IssmDouble* pviscosity, IssmDouble* pepsilon)=0; virtual void GetViscosity3d(IssmDouble* pviscosity3d, IssmDouble* pepsilon)=0; virtual void GetViscosity3dFS(IssmDouble* pviscosity3d, IssmDouble* epsilon)=0; Index: /issm/trunk-jpl/src/c/classes/Materials/Matice.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Materials/Matice.cpp (revision 17241) +++ /issm/trunk-jpl/src/c/classes/Materials/Matice.cpp (revision 17242) @@ -351,4 +351,69 @@ /*Build viscosity: viscosity=B/(2*A^e) */ A=exx*exx+eyy*eyy+2.*exy*exy; + if(A==0.){ + /*Maxiviscositym viscosity for 0 shear areas: */ + viscosity=2.5*2.e+17; + } + else{ + e=(n-1.)/(2.*n); + viscosity=(1.-D)*B/(2.*pow(A,e)); + } + } + } + + /*Checks in debugging mode*/ + if(viscosity<=0) _error_("Negative viscosity"); + _assert_(B>0); + _assert_(n>0); + _assert_(D>=0 && D<1); + + /*Return: */ + *pviscosity=viscosity; +} +/*}}}*/ +/*FUNCTION Matice::GetViscosity2dverticalHO {{{*/ +void Matice::GetViscosity2dverticalHO(IssmDouble* pviscosity, IssmDouble* epsilon){ + /*From a string tensor and a material object, return viscosity, using Glen's flow law. + (1-D) B + viscosity= -------------------------------------- + 2[ 2exx^2 + 2exy^2]^[(n-1)/2n] + + where viscosity is the viscotiy, B the flow law parameter , (u,v) the velocity + vector, and n the flow law exponent. + + If epsilon is NULL, it means this is the first time SystemMatrices is being run, and we + return 10^14, initial viscosity. + */ + + /*output: */ + IssmDouble viscosity; + + /*input strain rate: */ + IssmDouble exx,exy; + + /*Intermediary: */ + IssmDouble A,e; + IssmDouble B,D,n; + + /*Get B and n*/ + B=GetB(); + n=GetN(); + D=GetD(); + + if (n==1){ + /*Viscous behaviour! viscosity=B: */ + viscosity=(1-D)*B/2; + } + else{ + if((epsilon[0]==0) && (epsilon[1]==0) && (epsilon[2]==0)){ + viscosity=0.5*pow(10.,14); + } + else{ + /*Retrive strain rate components: */ + exx=epsilon[0]; + exy=epsilon[1]; + + /*Build viscosity: viscosity=B/(2*A^e) */ + A=2*exx*exx+2.*exy*exy; if(A==0.){ /*Maxiviscositym viscosity for 0 shear areas: */ Index: /issm/trunk-jpl/src/c/classes/Materials/Matice.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Materials/Matice.h (revision 17241) +++ /issm/trunk-jpl/src/c/classes/Materials/Matice.h (revision 17242) @@ -54,4 +54,5 @@ void GetViscosity2d(IssmDouble* pviscosity, IssmDouble* pepsilon); void GetViscosity2dvertical(IssmDouble* pviscosity, IssmDouble* pepsilon); + void GetViscosity2dverticalHO(IssmDouble* pviscosity, IssmDouble* pepsilon); void GetViscosity3d(IssmDouble* pviscosity3d, IssmDouble* pepsilon); void GetViscosity3dFS(IssmDouble* pviscosity3d, IssmDouble* epsilon); Index: /issm/trunk-jpl/src/c/classes/Materials/Matpar.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Materials/Matpar.h (revision 17241) +++ /issm/trunk-jpl/src/c/classes/Materials/Matpar.h (revision 17242) @@ -84,4 +84,5 @@ void GetViscosity2d(IssmDouble* pviscosity, IssmDouble* pepsilon){_error_("not supported");}; void GetViscosity2dvertical(IssmDouble* pviscosity, IssmDouble* pepsilon){_error_("not supported");}; + void GetViscosity2dverticalHO(IssmDouble* pviscosity, IssmDouble* pepsilon){_error_("not supported");}; void GetViscosity3d(IssmDouble* pviscosity3d, IssmDouble* pepsilon){_error_("not supported");}; void GetViscosity3dFS(IssmDouble* pviscosity3d, IssmDouble* epsilon){_error_("not supported");};