Index: /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp (revision 18292) +++ /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp (revision 18293) @@ -8,4 +8,5 @@ //#define FSANALYTICAL 12 +//#define LATERALFRICTION 1 /*Model processing*/ @@ -223,4 +224,7 @@ iomodel->FetchDataToInput(elements,LoadingforceXEnum); iomodel->FetchDataToInput(elements,LoadingforceYEnum); + #ifdef LATERALFRICTION + iomodel->FetchDataToInput(elements,MeshVertexonboundaryEnum); + #endif if(isdamage)iomodel->FetchDataToInput(elements,DamageDEnum); @@ -1218,5 +1222,11 @@ ElementMatrix* Ke1=CreateKMatrixSSAViscous(basalelement); ElementMatrix* Ke2=CreateKMatrixSSAFriction(basalelement); + #ifdef LATERALFRICTION + ElementMatrix* Ke3=CreateKMatrixSSALateralFriction(basalelement); + ElementMatrix* Ke =new ElementMatrix(Ke1,Ke2,Ke3); + delete Ke3; + #else ElementMatrix* Ke =new ElementMatrix(Ke1,Ke2); + #endif /*clean-up and return*/ @@ -1306,4 +1316,70 @@ delete friction; xDelete(xyz_list); + xDelete(B); + xDelete(D); + return Ke; +}/*}}}*/ +ElementMatrix* StressbalanceAnalysis::CreateKMatrixSSALateralFriction(Element* element){/*{{{*/ + + /*Return if element is inactive*/ + if(!element->IsIceInElement()) return NULL; + + /*If no boundary, return NULL*/ + if(!element->IsFaceOnBoundary()) return NULL; + + /*Intermediaries*/ + IssmDouble alpha2; + IssmDouble Jdet; + int domaintype; + IssmDouble icefront; + IssmDouble *xyz_list = NULL; + IssmDouble *xyz_list_boundary = NULL; + + /*Get problem dimension*/ + element->FindParam(&domaintype,DomainTypeEnum); + if(domaintype==Domain2DverticalEnum) return NULL; //not supported yet + + /*Fetch number of nodes and dof for this finite element*/ + int dim = 2; + int numnodes = element->GetNumberOfNodes(); + int numdof = numnodes*dim; + + /*Initialize Element matrix and vectors*/ + ElementMatrix* Ke = element->NewElementMatrix(SSAApproximationEnum); + IssmDouble* B = xNew(dim*numdof); + IssmDouble* D = xNewZeroInit(dim*dim); + + /*Retrieve all inputs and parameters*/ + element->GetVerticesCoordinates(&xyz_list); + element->GetLevelCoordinates(&xyz_list_boundary,xyz_list,MeshVertexonboundaryEnum,1.); + Input* icelevelset_input = element->GetInput(MaskIceLevelsetEnum); _assert_(icelevelset_input); + + /* Start looping on the number of gaussian points: */ + Gauss* gauss=element->NewGauss(xyz_list,xyz_list_boundary,3); + for(int ig=gauss->begin();igend();ig++){ + gauss->GaussPoint(ig); + + this->GetBSSAFriction(B,element,dim,xyz_list,gauss); + element->JacobianDeterminantSurface(&Jdet,xyz_list_boundary,gauss); + icelevelset_input->GetInputValue(&icefront, gauss); + if(icefront==0.) + alpha2=0.; + else + alpha2=2.e+12; + for(int i=0;iweight*Jdet; + + TripleMultiply(B,dim,numdof,1, + D,dim,dim,0, + B,dim,numdof,0, + &Ke->values[0],1); + } + + /*Transform Coordinate System*/ + if(dim==2) element->TransformStiffnessMatrixCoord(Ke,XYEnum); + + /*Clean up and return*/ + delete gauss; + xDelete(xyz_list); + xDelete(xyz_list_boundary); xDelete(B); xDelete(D); Index: /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.h =================================================================== --- /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.h (revision 18292) +++ /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.h (revision 18293) @@ -37,4 +37,5 @@ ElementMatrix* CreateKMatrixSSAViscous(Element* element); ElementMatrix* CreateKMatrixSSAFriction(Element* element); + ElementMatrix* CreateKMatrixSSALateralFriction(Element* element); ElementVector* CreatePVectorSSA(Element* element); ElementVector* CreatePVectorSSADrivingStress(Element* element); Index: /issm/trunk-jpl/src/c/classes/Elements/Element.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Element.h (revision 18292) +++ /issm/trunk-jpl/src/c/classes/Elements/Element.h (revision 18293) @@ -256,4 +256,5 @@ virtual void ZeroLevelsetCoordinates(IssmDouble** pxyz_zero,IssmDouble* xyz_list,int levelsetenum)=0; virtual void GetIcefrontCoordinates(IssmDouble** pxyz_front,IssmDouble* xyz_list,int levelsetenum)=0; + virtual void GetLevelCoordinates(IssmDouble** pxyz_front,IssmDouble* xyz_list,int levelsetenum,IssmDouble level)=0; virtual void AverageOntoPartition(Vector* partition_contributions,Vector* partition_areas,IssmDouble* vertex_response,IssmDouble* qmu_part)=0; Index: /issm/trunk-jpl/src/c/classes/Elements/Penta.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Penta.h (revision 18292) +++ /issm/trunk-jpl/src/c/classes/Elements/Penta.h (revision 18293) @@ -92,4 +92,5 @@ void ZeroLevelsetCoordinates(IssmDouble** pxyz_zero,IssmDouble* xyz_list,int levelsetenum); void GetIcefrontCoordinates(IssmDouble** pxyz_front,IssmDouble* xyz_list,int levelsetenum){_error_("not implemented yet");}; + void GetLevelCoordinates(IssmDouble** pxyz_front,IssmDouble* xyz_list,int levelsetenum,IssmDouble level){_error_("not implemented yet");}; void PositiveDegreeDay(IssmDouble* pdds,IssmDouble* pds,IssmDouble signorm); void ReduceMatrices(ElementMatrix* Ke,ElementVector* pe); Index: /issm/trunk-jpl/src/c/classes/Elements/Seg.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Seg.h (revision 18292) +++ /issm/trunk-jpl/src/c/classes/Elements/Seg.h (revision 18293) @@ -123,7 +123,8 @@ bool IsZeroLevelset(int levelset_enum){_error_("not implemented");}; bool IsIcefront(void); - bool IsFaceOnBoundary(void){_error_("not implemented yet");}; + bool IsFaceOnBoundary(void){_error_("not implemented yet");}; void ZeroLevelsetCoordinates(IssmDouble** pxyz_zero,IssmDouble* xyz_list,int levelsetenum){_error_("not implemented");}; void GetIcefrontCoordinates(IssmDouble** pxyz_front,IssmDouble* xyz_list,int levelsetenum); + void GetLevelCoordinates(IssmDouble** pxyz_front,IssmDouble* xyz_list,int levelsetenum,IssmDouble level){_error_("not implemented");}; void GetSolutionFromInputsOneDof(Vector* solution,int enum_type){_error_("not implemented yet");}; Index: /issm/trunk-jpl/src/c/classes/Elements/Tetra.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Tetra.h (revision 18292) +++ /issm/trunk-jpl/src/c/classes/Elements/Tetra.h (revision 18293) @@ -128,7 +128,8 @@ bool IsZeroLevelset(int levelset_enum){_error_("not implemented");}; bool IsIcefront(void); - bool IsFaceOnBoundary(void){_error_("not implemented yet");}; + bool IsFaceOnBoundary(void){_error_("not implemented yet");}; void ZeroLevelsetCoordinates(IssmDouble** pxyz_zero,IssmDouble* xyz_list,int levelsetenum); void GetIcefrontCoordinates(IssmDouble** pxyz_front,IssmDouble* xyz_list,int levelsetenum){_error_("not implemented yet");}; + void GetLevelCoordinates(IssmDouble** pxyz_front,IssmDouble* xyz_list,int levelsetenum,IssmDouble level){_error_("not implemented yet");}; void GetSolutionFromInputsOneDof(Vector* solution,int enum_type){_error_("not implemented yet");}; void InputDepthAverageAtBase(int enum_type,int average_enum_type){_error_("not implemented yet");}; Index: /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp (revision 18292) +++ /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp (revision 18293) @@ -846,4 +846,44 @@ for(i=0;i=0.){ + indicesfront[nrfrontnodes]=i; + nrfrontnodes++; + } + } + + _assert_(nrfrontnodes==2); + + /* arrange order of frontnodes such that they are oriented counterclockwise */ + if((NUMVERTICES+indicesfront[0]-indicesfront[1])%NUMVERTICES!=NUMVERTICES-1){ + int index=indicesfront[0]; + indicesfront[0]=indicesfront[1]; + indicesfront[1]=index; + } + + IssmDouble* xyz_front = xNew(3*nrfrontnodes); + /* Return nodes */ + for(i=0;i(indicesfront); +}/*}}}*/ +void Tria::GetLevelCoordinates(IssmDouble** pxyz_front,IssmDouble* xyz_list,int levelsetenum,IssmDouble level){/*{{{*/ + + /* Intermediaries */ + int i, dir,nrfrontnodes; + IssmDouble levelset[NUMVERTICES]; + + /*Recover parameters and values*/ + GetInputListOnVertices(&levelset[0],levelsetenum); + + int* indicesfront = xNew(NUMVERTICES); + /* Get nodes where there is no ice */ + nrfrontnodes=0; + for(i=0;i* partition_contributions,Vector* partition_areas,IssmDouble* vertex_response,IssmDouble* qmu_part);