Index: ../trunk-jpl/src/c/classes/Materials/Matice.cpp =================================================================== --- ../trunk-jpl/src/c/classes/Materials/Matice.cpp (revision 16349) +++ ../trunk-jpl/src/c/classes/Materials/Matice.cpp (revision 16350) @@ -291,6 +291,72 @@ *pviscosity=viscosity; } /*}}}*/ +/*FUNCTION Matice::GetViscosity2dvertical {{{*/ +void Matice::GetViscosity2dvertical(IssmDouble* pviscosity, IssmDouble* epsilon){ + /*From a string tensor and a material object, return viscosity, using Glen's flow law. + (1-D) B + viscosity= -------------------------------------- + 2[ exx^2+eyy^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,eyy,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]; + eyy=epsilon[1]; + exy=epsilon[2]; + + /*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::GetViscosity3d {{{*/ void Matice::GetViscosity3d(IssmDouble* pviscosity3d, IssmDouble* epsilon){ @@ -792,10 +858,10 @@ if(iomodel->meshtype==Mesh2DhorizontalEnum){ /*Intermediaries*/ - const int num_vertices = 3; //Tria has 3 vertices - IssmDouble nodeinputs[num_vertices]; - IssmDouble cmmininputs[num_vertices]; - IssmDouble cmmaxinputs[num_vertices]; + const int num_vertices = 3; //Tria has 3 vertices + IssmDouble nodeinputs[num_vertices]; + IssmDouble cmmininputs[num_vertices]; + IssmDouble cmmaxinputs[num_vertices]; /*Get B*/ if (iomodel->Data(MaterialsRheologyBEnum)) { @@ -866,7 +932,7 @@ /*Get D:*/ if (iomodel->Data(DamageDEnum)) { for(i=0;iData(DamageDEnum)[iomodel->elements[num_vertices*index+i]-1]; - this->inputs->AddInput(new TriaInput(DamageDbarEnum,nodeinputs,P1Enum)); + this->inputs->AddInput(new TriaInput(DamageDEnum,nodeinputs,P1Enum)); } } #ifdef _HAVE_3D_ Index: ../trunk-jpl/src/c/classes/Materials/Material.h =================================================================== --- ../trunk-jpl/src/c/classes/Materials/Material.h (revision 16349) +++ ../trunk-jpl/src/c/classes/Materials/Material.h (revision 16350) @@ -26,6 +26,7 @@ virtual void Configure(Elements* elements)=0; virtual void GetVectorFromInputs(Vector* vector,int input_enum)=0; virtual void GetViscosity2d(IssmDouble* pviscosity, IssmDouble* pepsilon)=0; + virtual void GetViscosity2dvertical(IssmDouble* pviscosity, IssmDouble* pepsilon)=0; virtual void GetViscosity3d(IssmDouble* pviscosity3d, IssmDouble* pepsilon)=0; virtual void GetViscosity3dFS(IssmDouble* pviscosity3d, IssmDouble* epsilon)=0; virtual void GetViscosityComplement(IssmDouble* pviscosity_complement, IssmDouble* pepsilon)=0; Index: ../trunk-jpl/src/c/classes/Materials/Matice.h =================================================================== --- ../trunk-jpl/src/c/classes/Materials/Matice.h (revision 16349) +++ ../trunk-jpl/src/c/classes/Materials/Matice.h (revision 16350) @@ -53,6 +53,7 @@ void GetVectorFromInputs(Vector* vector,int input_enum); void SetCurrentConfiguration(Elements* elements,Loads* loads,Nodes* nodes,Vertices* vertices,Materials* materials,Parameters* parameters); void GetViscosity2d(IssmDouble* pviscosity, IssmDouble* pepsilon); + void GetViscosity2dvertical(IssmDouble* pviscosity, IssmDouble* pepsilon); void GetViscosity3d(IssmDouble* pviscosity3d, IssmDouble* pepsilon); void GetViscosity3dFS(IssmDouble* pviscosity3d, IssmDouble* epsilon); void GetViscosityComplement(IssmDouble* pviscosity_complement, IssmDouble* pepsilon); Index: ../trunk-jpl/src/c/classes/Materials/Matpar.h =================================================================== --- ../trunk-jpl/src/c/classes/Materials/Matpar.h (revision 16349) +++ ../trunk-jpl/src/c/classes/Materials/Matpar.h (revision 16350) @@ -79,10 +79,11 @@ void InputUpdateFromIoModel(int index, IoModel* iomodel){_error_("not implemented yet");}; /*}}}*/ /*Material virtual functions resolution: {{{*/ - void InputDuplicate(int original_enum,int new_enum); - void Configure(Elements* elements); - void GetVectorFromInputs(Vector* vector,int input_enum){return;} + void InputDuplicate(int original_enum,int new_enum); + void Configure(Elements* elements); + void GetVectorFromInputs(Vector* vector,int input_enum){return;} void GetViscosity2d(IssmDouble* pviscosity, IssmDouble* pepsilon){_error_("not supported");}; + void GetViscosity2dvertical(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");}; void GetViscosityComplement(IssmDouble* pviscosity_complement, IssmDouble* pepsilon){_error_("not supported");}; Index: ../trunk-jpl/src/c/classes/Elements/Tria.cpp =================================================================== --- ../trunk-jpl/src/c/classes/Elements/Tria.cpp (revision 16349) +++ ../trunk-jpl/src/c/classes/Elements/Tria.cpp (revision 16350) @@ -19,7 +19,9 @@ /*}}}*/ /*Element macros*/ -#define NUMVERTICES 3 +#define NUMVERTICES 3 +#define NUMVERTICES1D 2 + /*Constructors/destructor/copy*/ /*FUNCTION Tria::Tria(){{{*/ Tria::Tria(){ @@ -218,7 +220,18 @@ switch(analysis_type){ #ifdef _HAVE_STRESSBALANCE_ case StressbalanceAnalysisEnum: - return CreateKMatrixStressbalanceSSA(); + int approximation; + inputs->GetInputValue(&approximation,ApproximationEnum); + switch(approximation){ + case SSAApproximationEnum: + return CreateKMatrixStressbalanceSSA(); + case FSApproximationEnum: + return CreateKMatrixStressbalanceFS(); + case NoneApproximationEnum: + return NULL; + default: + _error_("Approximation " << EnumToStringx(approximation) << " not supported yet"); + } break; case StressbalanceSIAAnalysisEnum: return CreateKMatrixStressbalanceSIA(); @@ -372,7 +385,18 @@ switch(analysis_type){ #ifdef _HAVE_STRESSBALANCE_ case StressbalanceAnalysisEnum: - return CreatePVectorStressbalanceSSA(); + int approximation; + inputs->GetInputValue(&approximation,ApproximationEnum); + switch(approximation){ + case SSAApproximationEnum: + return CreatePVectorStressbalanceSSA(); + case FSApproximationEnum: + return CreatePVectorStressbalanceFS(); + case NoneApproximationEnum: + return NULL; + default: + _error_("Approximation " << EnumToStringx(approximation) << " not supported yet"); + } break; case StressbalanceSIAAnalysisEnum: return CreatePVectorStressbalanceSIA(); @@ -1678,7 +1702,17 @@ switch(analysis_type){ #ifdef _HAVE_STRESSBALANCE_ case StressbalanceAnalysisEnum: - InputUpdateFromSolutionStressbalanceHoriz(solution); + int approximation; + inputs->GetInputValue(&approximation,ApproximationEnum); + if(approximation==FSApproximationEnum || approximation==NoneApproximationEnum){ + InputUpdateFromSolutionStressbalanceFS(solution); + } + else if (approximation==SSAApproximationEnum || approximation==SIAApproximationEnum){ + InputUpdateFromSolutionStressbalanceHoriz(solution); + } + else{ + _error_("approximation not supported yet"); + } break; case StressbalanceSIAAnalysisEnum: InputUpdateFromSolutionStressbalanceHoriz(solution); @@ -2055,6 +2089,50 @@ return onbed; } /*}}}*/ +/*FUNCTION Tria::HasEdgeOnBed {{{*/ +bool Tria::HasEdgeOnBed(){ + + IssmDouble values[NUMVERTICES]; + IssmDouble sum; + + /*Retrieve all inputs and parameters*/ + GetInputListOnVertices(&values[0],MeshVertexonbedEnum); + sum = values[0]+values[1]+values[2]; + + _assert_(sum==0. || sum==1. || sum==2.); + + if(sum>1.){ + return true; + } + else{ + return false; + } +} +/*}}}*/ +/*FUNCTION Tria::EdgeOnBedIndices{{{*/ +void Tria::EdgeOnBedIndices(int* pindex1,int* pindex2){ + + bool found=false; + IssmDouble values[NUMVERTICES]; + + /*Retrieve all inputs and parameters*/ + GetInputListOnVertices(&values[0],MeshVertexonbedEnum); + + for(int i=0;iNumberofNodesVelocity(); + int pnumnodes = this->NumberofNodesPressure(); + int numdof = vnumnodes*NDOF2 + pnumnodes*NDOF1; + + /*Prepare coordinate system list*/ + int* cs_list = xNew(vnumnodes+pnumnodes); + for(i=0;iparameters,FSvelocityEnum); + IssmDouble* B = xNew(5*numdof); + IssmDouble* Bprime = xNew(5*numdof); + IssmDouble* D = xNewZeroInit(5*5); + + /*Retrieve all inputs and parameters*/ + GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES); + parameters->FindParam(&FSreconditioning,StressbalanceFSreconditioningEnum); + Input* vx_input=inputs->GetInput(VxEnum); _assert_(vx_input); + Input* vy_input=inputs->GetInput(VyEnum); _assert_(vy_input); + + /* Start looping on the number of gaussian points: */ + gauss=new GaussTria(5); + for(int ig=gauss->begin();igend();ig++){ + + gauss->GaussPoint(ig); + + GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss); + GetBFS(B,&xyz_list[0][0],gauss); + GetBprimeFS(Bprime,&xyz_list[0][0],gauss); + + this->GetStrainRate2d(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input); + material->GetViscosity2dvertical(&viscosity,&epsilon[0]); + + D_scalar=gauss->weight*Jdet; + for(i=0;i<3;i++) D[i*5+i] = +D_scalar*2.*viscosity; + for(i=3;i<5;i++) D[i*5+i] = -D_scalar*FSreconditioning; + + TripleMultiply(B,5,numdof,1, + D,5,5,0, + Bprime,5,numdof,0, + Ke->values,1); + } + + /*Transform Coordinate System*/ + TransformStiffnessMatrixCoord(Ke,nodes,(vnumnodes+pnumnodes),cs_list); + + /*Clean up and return*/ + delete gauss; + xDelete(cs_list); + xDelete(B); + xDelete(Bprime); + xDelete(D); + return Ke; +} +/*}}}*/ +/*FUNCTION Tria::CreateKMatrixStressbalanceFSFriction{{{*/ +ElementMatrix* Tria::CreateKMatrixStressbalanceFSFriction(void){ + + /*Intermediaries */ + int i,j; + int analysis_type,approximation; + IssmDouble alpha2,Jdet2d; + IssmDouble FSreconditioning,viscosity; + IssmDouble epsilon[3]; /* epsilon=[exx,eyy,exy];*/ + IssmDouble xyz_list[NUMVERTICES][3]; + IssmDouble xyz_list_seg[NUMVERTICES1D][3]; + Friction *friction = NULL; + GaussTria *gauss = NULL; + + /*Initialize Element matrix and return if necessary*/ + if(IsFloating() || !HasEdgeOnBed()) return NULL; + + _error_("STOP"); +} +/*}}}*/ /*FUNCTION Tria::CreateKMatrixStressbalanceSSA {{{*/ ElementMatrix* Tria::CreateKMatrixStressbalanceSSA(void){ @@ -3225,6 +3409,183 @@ return Ke; } /*}}}*/ +/*FUNCTION Tria::CreatePVectorStressbalanceFS {{{*/ +ElementVector* Tria::CreatePVectorStressbalanceFS(void){ + + ElementVector* pe1; + ElementVector* pe2; + ElementVector* pe3; + ElementVector* pe; + + /*compute all stiffness matrices for this element*/ + pe1=CreatePVectorStressbalanceFSViscous(); + pe2=CreatePVectorStressbalanceFSShelf(); + pe3=CreatePVectorStressbalanceFSFront(); + pe =new ElementVector(pe1,pe2,pe3); + + /*clean-up and return*/ + delete pe1; + delete pe2; + delete pe3; + return pe; +} +/*}}}*/ +/*FUNCTION Tria::CreatePVectorStressbalanceFSFront{{{*/ +ElementVector* Tria::CreatePVectorStressbalanceFSFront(void){ + + /*Intermediaries */ + int i; + IssmDouble ls[NUMVERTICES]; + IssmDouble xyz_list[NUMVERTICES][3]; + bool isfront; + + /*Retrieve all inputs and parameters*/ + GetInputListOnVertices(&ls[0],MaskIceLevelsetEnum); + + /*If the level set is awlays <=0, there is no ice front here*/ + isfront = false; + if(ls[0]>0. || ls[1]>0. || ls[2]>0.){ + if(ls[0]*ls[1]<0. || ls[0]*ls[2]<0. || (ls[0]*ls[1]+ls[0]*ls[2]+ls[1]*ls[2]==0.)){ + isfront = true; + } + } + + /*If no front, return NULL*/ + if(!isfront) return NULL; + + _error_("STOP"); +} +/*}}}*/ +/*FUNCTION Tria::CreatePVectorStressbalanceFSViscous {{{*/ +ElementVector* Tria::CreatePVectorStressbalanceFSViscous(void){ + + /*Intermediaries*/ + int i; + int approximation; + IssmDouble Jdet,gravity,rho_ice; + IssmDouble forcex,forcey; + IssmDouble xyz_list[NUMVERTICES][3]; + GaussTria *gauss=NULL; + + /*Fetch number of nodes and dof for this finite element*/ + int vnumnodes = this->NumberofNodesVelocity(); + int pnumnodes = this->NumberofNodesPressure(); + + /*Prepare coordinate system list*/ + int* cs_list = xNew(vnumnodes+pnumnodes); + for(i=0;iparameters,FSvelocityEnum); + IssmDouble* vbasis = xNew(vnumnodes); + + /*Retrieve all inputs and parameters*/ + GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES); + Input* loadingforcex_input=inputs->GetInput(LoadingforceXEnum); _assert_(loadingforcex_input); + Input* loadingforcey_input=inputs->GetInput(LoadingforceYEnum); _assert_(loadingforcey_input); + rho_ice = matpar->GetRhoIce(); + gravity = matpar->GetG(); + + /* Start looping on the number of gaussian points: */ + gauss=new GaussTria(5); + for(int ig=gauss->begin();igend();ig++){ + + gauss->GaussPoint(ig); + + GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss); + GetNodalFunctionsVelocity(vbasis, gauss); + + loadingforcex_input->GetInputValue(&forcex,gauss); + loadingforcey_input->GetInputValue(&forcey,gauss); + + for(i=0;ivalues[i*NDOF2+1] += -rho_ice*gravity*Jdet*gauss->weight*vbasis[i]; + pe->values[i*NDOF2+0] += +rho_ice*forcex *Jdet*gauss->weight*vbasis[i]; + pe->values[i*NDOF2+1] += +rho_ice*forcey *Jdet*gauss->weight*vbasis[i]; + } + } + + /*Transform coordinate system*/ + TransformLoadVectorCoord(pe,nodes,vnumnodes+pnumnodes,cs_list); + + /*Clean up and return*/ + delete gauss; + xDelete(cs_list); + xDelete(vbasis); + return pe; +} +/*}}}*/ +/*FUNCTION Tria::CreatePVectorStressbalanceFSShelf{{{*/ +ElementVector* Tria::CreatePVectorStressbalanceFSShelf(void){ + + /*Intermediaries*/ + int i,j; + int indices[2]; + IssmDouble gravity,rho_water,bed,water_pressure; + IssmDouble normal_vel,vx,vy,dt; + IssmDouble xyz_list_seg[NUMVERTICES1D][3]; + IssmDouble xyz_list[NUMVERTICES][3]; + IssmDouble normal[2]; + IssmDouble Jdet; + + /*Initialize Element vector and return if necessary*/ + if(!HasEdgeOnBed() || !IsFloating()) return NULL; + + /*Fetch number of nodes and dof for this finite element*/ + int vnumnodes = this->NumberofNodesVelocity(); + int pnumnodes = this->NumberofNodesPressure(); + + /*Prepare coordinate system list*/ + int* cs_list = xNew(vnumnodes+pnumnodes); + for(i=0;iparameters,FSvelocityEnum); + IssmDouble* vbasis = xNew(vnumnodes); + + /*Retrieve all inputs and parameters*/ + rho_water=matpar->GetRhoWater(); + gravity=matpar->GetG(); + GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES); + Input* bed_input=inputs->GetInput(BedEnum); _assert_(bed_input); + + /*Get vertex indices that lie on bed*/ + this->EdgeOnBedIndices(&indices[0],&indices[1]); + + for(i=0;ibegin();igend();ig++){ + + gauss->GaussPoint(ig); + + GetSegmentJacobianDeterminant(&Jdet,&xyz_list_seg[0][0],gauss); + GetNodalFunctionsVelocity(vbasis, gauss); + + GetSegmentNormal(&normal[0],xyz_list_seg); + _assert_(normal[1]<0.); + bed_input->GetInputValue(&bed, gauss); + water_pressure=gravity*rho_water*bed; + + for(i=0;ivalues[2*i+0]+=(water_pressure)*gauss->weight*Jdet*vbasis[i]*normal[0]; + pe->values[2*i+1]+=(water_pressure)*gauss->weight*Jdet*vbasis[i]*normal[1]; + } + } + + /*Transform coordinate system*/ + TransformLoadVectorCoord(pe,nodes,(vnumnodes+pnumnodes),cs_list); + + /*Clean up and return*/ + delete gauss; + xDelete(cs_list); + xDelete(vbasis); + return pe; +} +/*}}}*/ /*FUNCTION Tria::CreatePVectorStressbalanceSSA {{{*/ ElementVector* Tria::CreatePVectorStressbalanceSSA(){ @@ -3695,6 +4056,79 @@ } /*}}}*/ +/*FUNCTION Tria::InputUpdateFromSolutionStressbalanceFS {{{*/ +void Tria::InputUpdateFromSolutionStressbalanceFS(IssmDouble* solution){ + + int i; + int* vdoflist=NULL; + int* pdoflist=NULL; + IssmDouble FSreconditioning; + + /*Fetch number of nodes and dof for this finite element*/ + int vnumnodes = this->NumberofNodesVelocity(); + int pnumnodes = this->NumberofNodesPressure(); + int vnumdof = vnumnodes*NDOF2; + int pnumdof = pnumnodes*NDOF1; + + /*Initialize values*/ + IssmDouble* vvalues = xNew(vnumdof); + IssmDouble* pvalues = xNew(pnumdof); + IssmDouble* vx = xNew(vnumnodes); + IssmDouble* vy = xNew(vnumnodes); + IssmDouble* vel = xNew(vnumnodes); + IssmDouble* pressure = xNew(pnumnodes); + + /*Get dof list: */ + GetDofListVelocity(&vdoflist,GsetEnum); + GetDofListPressure(&pdoflist,GsetEnum); + + /*Use the dof list to index into the solution vector: */ + for(i=0;i(vx[i])) _error_("NaN found in solution vector"); + if(xIsNan(vy[i])) _error_("NaN found in solution vector"); + } + for(i=0;i(pressure[i])) _error_("NaN found in solution vector"); + } + + /*Recondition pressure and compute vel: */ + this->parameters->FindParam(&FSreconditioning,StressbalanceFSreconditioningEnum); + for(i=0;iinputs->ChangeEnum(VxEnum,VxPicardEnum); + this->inputs->ChangeEnum(VyEnum,VyPicardEnum); + this->inputs->ChangeEnum(PressureEnum,PressurePicardEnum); + + /*Add vx and vy as inputs to the tria element: */ + this->inputs->AddInput(new TriaInput(VxEnum,vx,P1Enum)); + this->inputs->AddInput(new TriaInput(VyEnum,vy,P1Enum)); + this->inputs->AddInput(new TriaInput(VelEnum,vel,P1Enum)); + this->inputs->AddInput(new TriaInput(PressureEnum,pressure,P1Enum)); + + /*Free ressources:*/ + xDelete(pressure); + xDelete(vel); + xDelete(vy); + xDelete(vx); + xDelete(vvalues); + xDelete(pvalues); + xDelete(vdoflist); + xDelete(pdoflist); +} +/*}}}*/ /*FUNCTION Tria::InputUpdateFromSolutionStressbalanceSIA {{{*/ void Tria::InputUpdateFromSolutionStressbalanceSIA(IssmDouble* solution){ Index: ../trunk-jpl/src/c/classes/Elements/TriaRef.h =================================================================== --- ../trunk-jpl/src/c/classes/Elements/TriaRef.h (revision 16349) +++ ../trunk-jpl/src/c/classes/Elements/TriaRef.h (revision 16350) @@ -22,8 +22,10 @@ void SetElementType(int type,int type_counter); /*Numerics*/ + void GetBFS(IssmDouble* B_prime, IssmDouble* xyz_list, GaussTria* gauss); void GetBSSA(IssmDouble* B, IssmDouble* xyz_list, GaussTria* gauss); void GetBSSAFS(IssmDouble* B , IssmDouble* xyz_list, GaussTria* gauss); + void GetBprimeFS(IssmDouble* B_prime, IssmDouble* xyz_list, GaussTria* gauss); void GetBprimeSSA(IssmDouble* Bprime, IssmDouble* xyz_list, GaussTria* gauss); void GetBprimeSSAFS(IssmDouble* Bprime, IssmDouble* xyz_list, GaussTria* gauss); void GetBprimeMasstransport(IssmDouble* Bprime, IssmDouble* xyz_list, GaussTria* gauss); @@ -35,10 +37,14 @@ void GetJacobianDeterminant(IssmDouble* Jdet, IssmDouble* xyz_list,GaussTria* gauss); void GetJacobianInvert(IssmDouble* Jinv, IssmDouble* xyz_list,GaussTria* gauss); void GetNodalFunctions(IssmDouble* basis,GaussTria* gauss); + void GetNodalFunctionsVelocity(IssmDouble* basis, GaussTria* gauss); + void GetNodalFunctionsPressure(IssmDouble* basis, GaussTria* gauss); void GetSegmentNodalFunctions(IssmDouble* basis,GaussTria* gauss, int index1,int index2); void GetSegmentBFlux(IssmDouble* B,GaussTria* gauss, int index1,int index2); void GetSegmentBprimeFlux(IssmDouble* Bprime,GaussTria* gauss, int index1,int index2); void GetNodalFunctionsDerivatives(IssmDouble* dbasis,IssmDouble* xyz_list, GaussTria* gauss); + void GetNodalFunctionsDerivativesVelocity(IssmDouble* dbasis,IssmDouble* xyz_list,GaussTria* gauss); + void GetNodalFunctionsDerivativesPressure(IssmDouble* dbasis,IssmDouble* xyz_list,GaussTria* gauss); void GetNodalFunctionsDerivativesReference(IssmDouble* dbasis,GaussTria* gauss); void GetInputValue(IssmDouble* pp, IssmDouble* plist, GaussTria* gauss); void GetInputDerivativeValue(IssmDouble* pp, IssmDouble* plist,IssmDouble* xyz_list, GaussTria* gauss); Index: ../trunk-jpl/src/c/classes/Elements/PentaRef.cpp =================================================================== --- ../trunk-jpl/src/c/classes/Elements/PentaRef.cpp (revision 16349) +++ ../trunk-jpl/src/c/classes/Elements/PentaRef.cpp (revision 16350) @@ -333,9 +333,9 @@ /*Compute B matrix. B=[Bv1 Bv2 ... Bp1 Bp2 ...] where Bvi is of size 3*NDOF3. * For node i, Bvi can be expressed in the actual coordinate system - * by: Bvi=[ dh/dx 0 0 ] + * by: Bvi=[ dh/dx 0 0 ] * [ 0 dh/dy 0 ] - * [ 0 0 dh/dy ] + * [ 0 0 dh/dz ] * [ 1/2*dh/dy 1/2*dh/dx 0 ] * [ 1/2*dh/dz 0 1/2*dh/dx ] * [ 0 1/2*dh/dz 1/2*dh/dy ] Index: ../trunk-jpl/src/c/classes/Elements/Tria.h =================================================================== --- ../trunk-jpl/src/c/classes/Elements/Tria.h (revision 16349) +++ ../trunk-jpl/src/c/classes/Elements/Tria.h (revision 16350) @@ -83,6 +83,8 @@ int GetNumberOfNodes(void); int Sid(); bool IsOnBed(); + bool HasEdgeOnBed(); + void EdgeOnBedIndices(int* pindex1,int* pindex); bool IsFloating(); bool IsNodeOnShelfFromFlags(IssmDouble* flags); bool NoIceInElement(); @@ -232,15 +234,24 @@ ElementMatrix* CreateKMatrixStressbalanceSSAViscous(void); ElementMatrix* CreateKMatrixStressbalanceSSAFriction(void); ElementMatrix* CreateKMatrixStressbalanceSIA(void); + ElementMatrix* CreateKMatrixStressbalanceFS(void); + ElementMatrix* CreateKMatrixStressbalanceFSViscous(void); + ElementMatrix* CreateKMatrixStressbalanceFSFriction(void); ElementVector* CreatePVectorStressbalanceSSA(void); ElementVector* CreatePVectorStressbalanceSSADrivingStress(void); ElementVector* CreatePVectorStressbalanceSSAFront(void); ElementVector* CreatePVectorStressbalanceSIA(void); + ElementVector* CreatePVectorStressbalanceFS(void); + ElementVector* CreatePVectorStressbalanceFSFront(void); + ElementVector* CreatePVectorStressbalanceFSViscous(void); + void PVectorGLSstabilization(ElementVector* pe); + ElementVector* CreatePVectorStressbalanceFSShelf(void); ElementMatrix* CreateJacobianStressbalanceSSA(void); void GetSolutionFromInputsStressbalanceFS(Vector* solution); void GetSolutionFromInputsStressbalanceHoriz(Vector* solution); void GetSolutionFromInputsStressbalanceSIA(Vector* solution); void InputUpdateFromSolutionStressbalanceHoriz( IssmDouble* solution); + void InputUpdateFromSolutionStressbalanceFS( IssmDouble* solution); void InputUpdateFromSolutionStressbalanceSIA( IssmDouble* solution); #endif Index: ../trunk-jpl/src/c/classes/Elements/Penta.cpp =================================================================== --- ../trunk-jpl/src/c/classes/Elements/Penta.cpp (revision 16349) +++ ../trunk-jpl/src/c/classes/Elements/Penta.cpp (revision 16350) @@ -10657,7 +10657,7 @@ IssmDouble* vy = xNew(vnumnodes); IssmDouble* vz = xNew(vnumnodes); IssmDouble* vel = xNew(vnumnodes); - IssmDouble* pressure = xNew(vnumnodes); + IssmDouble* pressure = xNew(pnumnodes); /*Get dof list: */ GetDofListVelocity(&vdoflist,GsetEnum); Index: ../trunk-jpl/src/c/classes/Elements/TriaRef.cpp =================================================================== --- ../trunk-jpl/src/c/classes/Elements/TriaRef.cpp (revision 16349) +++ ../trunk-jpl/src/c/classes/Elements/TriaRef.cpp (revision 16350) @@ -199,6 +199,116 @@ xDelete(basis); } /*}}}*/ +/*FUNCTION TriaRef::GetBFS {{{*/ +void TriaRef::GetBFS(IssmDouble* B, IssmDouble* xyz_list, GaussTria* gauss){ + /*Compute B matrix. B=[Bv1 Bv2 ... Bp1 Bp2 ...] where Bvi is of size 3*NDOF3. + * For node i, Bvi can be expressed in the actual coordinate system + * by: Bvi=[ dphi/dx 0 ] + * [ 0 dphi/dy ] + * [ 1/2*dphi/dy 1/2*dphi/dx] + * [ 0 0 ] + * [ dphi/dx dphi/dy ] + * + * by: Bpi=[ 0 ] + * [ 0 ] + * [ 0 ] + * [ phi_p ] + * [ 0 ] + * where phi is the interpolation function for node i. + * Same thing for Bb except the last column that does not exist. + */ + + /*Fetch number of nodes for this finite element*/ + int pnumnodes = this->NumberofNodesPressure(); + int vnumnodes = this->NumberofNodesVelocity(); + + /*Get nodal functions derivatives*/ + IssmDouble* vdbasis=xNew(2*vnumnodes); + IssmDouble* pbasis =xNew(pnumnodes); + GetNodalFunctionsDerivativesVelocity(vdbasis,xyz_list,gauss); + GetNodalFunctionsPressure(pbasis,gauss); + + /*Build B: */ + for(int i=0;i(vdbasis); + xDelete(pbasis); +} +/*}}}*/ +/*FUNCTION TriaRef::GetBprimeFS {{{*/ +void TriaRef::GetBprimeFS(IssmDouble* B_prime, IssmDouble* xyz_list, GaussTria* gauss){ + /* Compute B' matrix. B'=[B1' B2' B3' B4' B5' B6' Bb'] where Bi' is of size 3*NDOF2. + * For node i, Bi' can be expressed in the actual coordinate system + * by: + * Bvi' = [ dphi/dx 0 ] + * [ 0 dphi/dy ] + * [ dphi/dy dphi/dx ] + * [ dphi/dx dphi/dy ] + * [ 0 0 ] + * + * by: Bpi=[ 0 ] + * [ 0 ] + * [ 0 ] + * [ 0 ] + * [ phi ] + * where phi is the interpolation function for node i. + */ + + /*Fetch number of nodes for this finite element*/ + int pnumnodes = this->NumberofNodesPressure(); + int vnumnodes = this->NumberofNodesVelocity(); + + /*Get nodal functions derivatives*/ + IssmDouble* vdbasis=xNew(2*vnumnodes); + IssmDouble* pbasis =xNew(pnumnodes); + GetNodalFunctionsDerivativesVelocity(vdbasis,xyz_list,gauss); + GetNodalFunctionsPressure(pbasis,gauss); + + /*Build B_prime: */ + for(int i=0;i(vdbasis); + xDelete(pbasis); +} +/*}}}*/ /*FUNCTION TriaRef::GetBMasstransport{{{*/ void TriaRef::GetBMasstransport(IssmDouble* B, IssmDouble* xyz_list, GaussTria* gauss){ /*Compute B matrix. B=[B1 B2 B3] where Bi is of size 3*NDOF2. @@ -457,6 +567,76 @@ } } /*}}}*/ +/*FUNCTION TriaRef::GetNodalFunctionsVelocity{{{*/ +void TriaRef::GetNodalFunctionsVelocity(IssmDouble* basis,GaussTria* gauss){ + /*This routine returns the values of the nodal functions at the gaussian point.*/ + + switch(this->element_type){ + case P1P1Enum: + this->element_type = P1Enum; + this->GetNodalFunctions(basis,gauss); + this->element_type = P1P1Enum; + return; + case P1P1GLSEnum: + this->element_type = P1Enum; + this->GetNodalFunctions(basis,gauss); + this->element_type = P1P1GLSEnum; + return; + case MINIcondensedEnum: + this->element_type = P1bubbleEnum; + this->GetNodalFunctions(basis,gauss); + this->element_type = MINIcondensedEnum; + return; + case MINIEnum: + this->element_type = P1bubbleEnum; + this->GetNodalFunctions(basis,gauss); + this->element_type = MINIEnum; + return; + case TaylorHoodEnum: + this->element_type = P2Enum; + this->GetNodalFunctions(basis,gauss); + this->element_type = TaylorHoodEnum; + return; + default: + _error_("Element type "<element_type)<<" not supported yet"); + } +} +/*}}}*/ +/*FUNCTION TriaRef::GetNodalFunctionsPressure{{{*/ +void TriaRef::GetNodalFunctionsPressure(IssmDouble* basis,GaussTria* gauss){ + /*This routine returns the values of the nodal functions at the gaussian point.*/ + + switch(this->element_type){ + case P1P1Enum: + this->element_type = P1Enum; + this->GetNodalFunctions(basis,gauss); + this->element_type = P1P1Enum; + return; + case P1P1GLSEnum: + this->element_type = P1Enum; + this->GetNodalFunctions(basis,gauss); + this->element_type = P1P1GLSEnum; + return; + case MINIcondensedEnum: + this->element_type = P1Enum; + this->GetNodalFunctions(basis,gauss); + this->element_type = MINIcondensedEnum; + return; + case MINIEnum: + this->element_type = P1Enum; + this->GetNodalFunctions(basis,gauss); + this->element_type = MINIEnum; + return; + case TaylorHoodEnum: + this->element_type = P1Enum; + this->GetNodalFunctions(basis,gauss); + this->element_type = TaylorHoodEnum; + return; + default: + _error_("Element type "<element_type)<<" not supported yet"); + } +} +/*}}}*/ /*FUNCTION TriaRef::GetSegmentNodalFunctions{{{*/ void TriaRef::GetSegmentNodalFunctions(IssmDouble* basis,GaussTria* gauss,int index1,int index2){ /*This routine returns the values of the nodal functions at the gaussian point.*/ @@ -528,6 +708,72 @@ } /*}}}*/ +/*FUNCTION TriaRef::GetNodalFunctionsDerivativesPressure{{{*/ +void TriaRef::GetNodalFunctionsDerivativesPressure(IssmDouble* dbasis,IssmDouble* xyz_list, GaussTria* gauss){ + switch(this->element_type){ + case P1P1Enum: + this->element_type = P1Enum; + this->GetNodalFunctionsDerivatives(dbasis,xyz_list,gauss); + this->element_type = P1P1Enum; + return; + case P1P1GLSEnum: + this->element_type = P1Enum; + this->GetNodalFunctionsDerivatives(dbasis,xyz_list,gauss); + this->element_type = P1P1GLSEnum; + return; + case MINIcondensedEnum: + this->element_type = P1Enum; + this->GetNodalFunctionsDerivatives(dbasis,xyz_list,gauss); + this->element_type = MINIcondensedEnum; + return; + case MINIEnum: + this->element_type = P1Enum; + this->GetNodalFunctionsDerivatives(dbasis,xyz_list,gauss); + this->element_type = MINIEnum; + return; + case TaylorHoodEnum: + this->element_type = P1Enum; + this->GetNodalFunctionsDerivatives(dbasis,xyz_list,gauss); + this->element_type = TaylorHoodEnum; + return; + default: + _error_("Element type "<element_type)<<" not supported yet"); + } +} +/*}}}*/ +/*FUNCTION TriaRef::GetNodalFunctionsDerivativesVelocity{{{*/ +void TriaRef::GetNodalFunctionsDerivativesVelocity(IssmDouble* dbasis,IssmDouble* xyz_list, GaussTria* gauss){ + switch(this->element_type){ + case P1P1Enum: + this->element_type = P1Enum; + this->GetNodalFunctionsDerivatives(dbasis,xyz_list,gauss); + this->element_type = P1P1Enum; + return; + case P1P1GLSEnum: + this->element_type = P1Enum; + this->GetNodalFunctionsDerivatives(dbasis,xyz_list,gauss); + this->element_type = P1P1GLSEnum; + return; + case MINIcondensedEnum: + this->element_type = P1bubbleEnum; + this->GetNodalFunctionsDerivatives(dbasis,xyz_list,gauss); + this->element_type = MINIcondensedEnum; + return; + case MINIEnum: + this->element_type = P1bubbleEnum; + this->GetNodalFunctionsDerivatives(dbasis,xyz_list,gauss); + this->element_type = MINIEnum; + return; + case TaylorHoodEnum: + this->element_type = P2Enum; + this->GetNodalFunctionsDerivatives(dbasis,xyz_list,gauss); + this->element_type = TaylorHoodEnum; + return; + default: + _error_("Element type "<element_type)<<" not supported yet"); + } +} +/*}}}*/ /*FUNCTION TriaRef::GetNodalFunctionsDerivativesReference{{{*/ void TriaRef::GetNodalFunctionsDerivativesReference(IssmDouble* dbasis,GaussTria* gauss){ /*This routine returns the values of the nodal functions derivatives (with respect to the