Index: /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp (revision 25513) +++ /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp (revision 25514) @@ -929,4 +929,6 @@ parameters->AddObject(iomodel->CopyConstantObject("md.flowequation.isFS",FlowequationIsFSEnum)); parameters->AddObject(iomodel->CopyConstantObject("md.flowequation.fe_FS",FlowequationFeFSEnum)); + parameters->AddObject(iomodel->CopyConstantObject("md.flowequation.isNitscheBC",FlowequationIsNitscheEnum)); + parameters->AddObject(iomodel->CopyConstantObject("md.flowequation.FSNitscheGamma",FeFSNitscheGammaEnum)); parameters->AddObject(iomodel->CopyConstantObject("md.stressbalance.restol",StressbalanceRestolEnum)); parameters->AddObject(iomodel->CopyConstantObject("md.stressbalance.reltol",StressbalanceReltolEnum)); @@ -3433,5 +3435,8 @@ /*Get type of algorithm*/ int fe_FS; + bool isNitsche; + element->FindParam(&fe_FS,FlowequationFeFSEnum); + element->FindParam(&isNitsche,FlowequationIsNitscheEnum); /*compute all stiffness matrices for this element*/ @@ -3441,8 +3446,16 @@ else if(fe_FS==LATaylorHoodEnum || fe_FS==LACrouzeixRaviartEnum) Ke1=CreateKMatrixFSViscousLA(element); + else if(fe_FS==P1P1GLSEnum) + Ke1=CreateKMatrixFSViscousGLS(element); else Ke1=CreateKMatrixFSViscous(element); - ElementMatrix* Ke2=CreateKMatrixFSFriction(element); + ElementMatrix* Ke2; + if (isNitsche) { + Ke2 = CreateKMatrixFSFrictionNitsche(element); + } + else { + Ke2 = CreateKMatrixFSFriction(element); + } ElementMatrix* Ke3=CreateKMatrixFSShelf(element); ElementMatrix* Ke =new ElementMatrix(Ke1,Ke2,Ke3); @@ -3564,22 +3577,21 @@ element->material->ViscosityFS(&viscosity,dim,xyz_list,gauss,vx_input,vy_input,vz_input); - if(dim==3){ + if(dim==2 || dim==3){ /*Stress balance*/ for(int i=0;ivalues[(3*i+0)*numdof+3*j+0] += gauss->weight*Jdet*viscosity*(2.*vdbasis[0*vnumnodes+j]*vdbasis[0*vnumnodes+i] + vdbasis[1*vnumnodes+j]*vdbasis[1*vnumnodes+i] + vdbasis[2*vnumnodes+j]*vdbasis[2*vnumnodes+i]); - Ke->values[(3*i+0)*numdof+3*j+1] += gauss->weight*Jdet*viscosity*(vdbasis[0*vnumnodes+j]*vdbasis[1*vnumnodes+i]); - Ke->values[(3*i+0)*numdof+3*j+2] += gauss->weight*Jdet*viscosity*(vdbasis[0*vnumnodes+j]*vdbasis[2*vnumnodes+i]); - Ke->values[(3*i+1)*numdof+3*j+0] += gauss->weight*Jdet*viscosity*(vdbasis[1*vnumnodes+j]*vdbasis[0*vnumnodes+i]); - Ke->values[(3*i+1)*numdof+3*j+1] += gauss->weight*Jdet*viscosity*(vdbasis[0*vnumnodes+j]*vdbasis[0*vnumnodes+i] + 2.*vdbasis[1*vnumnodes+j]*vdbasis[1*vnumnodes+i] + vdbasis[2*vnumnodes+j]*vdbasis[2*vnumnodes+i]); - Ke->values[(3*i+1)*numdof+3*j+2] += gauss->weight*Jdet*viscosity*(vdbasis[1*vnumnodes+j]*vdbasis[2*vnumnodes+i]); - Ke->values[(3*i+2)*numdof+3*j+0] += gauss->weight*Jdet*viscosity*(vdbasis[2*vnumnodes+j]*vdbasis[0*vnumnodes+i]); - Ke->values[(3*i+2)*numdof+3*j+1] += gauss->weight*Jdet*viscosity*(vdbasis[2*vnumnodes+j]*vdbasis[1*vnumnodes+i]); - Ke->values[(3*i+2)*numdof+3*j+2] += gauss->weight*Jdet*viscosity*(vdbasis[0*vnumnodes+j]*vdbasis[0*vnumnodes+i] + vdbasis[1*vnumnodes+j]*vdbasis[1*vnumnodes+i] + 2.*vdbasis[2*vnumnodes+j]*vdbasis[2*vnumnodes+i]); + for (int p=0;pvalues[(dim*i+p)*numdof+dim*j+p] += gauss->weight*Jdet*viscosity*(vdbasis[q*vnumnodes+j]*vdbasis[q*vnumnodes+i]); + /* All the entries */ + Ke->values[(dim*i+p)*numdof+dim*j+q] += gauss->weight*Jdet*viscosity*(vdbasis[p*vnumnodes+j]*vdbasis[q*vnumnodes+i]); + } + } } for(int k=0;kvalues[(3*i+0)*numdof+3*vnumnodes+k] += gauss->weight*Jdet*FSreconditioning*(-pbasis[k]*vdbasis[0*vnumnodes+i]); - Ke->values[(3*i+1)*numdof+3*vnumnodes+k] += gauss->weight*Jdet*FSreconditioning*(-pbasis[k]*vdbasis[1*vnumnodes+i]); - Ke->values[(3*i+2)*numdof+3*vnumnodes+k] += gauss->weight*Jdet*FSreconditioning*(-pbasis[k]*vdbasis[2*vnumnodes+i]); + for (int p=0;pvalues[(dim*i+p)*numdof+dim*vnumnodes+k] += gauss->weight*Jdet*FSreconditioning*(-pbasis[k]*vdbasis[p*vnumnodes+i]); + } } } @@ -3587,29 +3599,7 @@ for(int k=0;kvalues[(3*vnumnodes+k)*numdof+3*j+0] += gauss->weight*Jdet*(-FSreconditioning*vdbasis[0*vnumnodes+j]*pbasis[k]); - Ke->values[(3*vnumnodes+k)*numdof+3*j+1] += gauss->weight*Jdet*(-FSreconditioning*vdbasis[1*vnumnodes+j]*pbasis[k]); - Ke->values[(3*vnumnodes+k)*numdof+3*j+2] += gauss->weight*Jdet*(-FSreconditioning*vdbasis[2*vnumnodes+j]*pbasis[k]); - } - } - } - else{ - /*Stress balance*/ - for(int i=0;ivalues[(2*i+0)*numdof+2*j+0] += gauss->weight*Jdet*viscosity*(2.*vdbasis[0*vnumnodes+j]*vdbasis[0*vnumnodes+i] + vdbasis[1*vnumnodes+j]*vdbasis[1*vnumnodes+i]); - Ke->values[(2*i+0)*numdof+2*j+1] += gauss->weight*Jdet*viscosity*(vdbasis[0*vnumnodes+j]*vdbasis[1*vnumnodes+i]); - Ke->values[(2*i+1)*numdof+2*j+0] += gauss->weight*Jdet*viscosity*(vdbasis[1*vnumnodes+j]*vdbasis[0*vnumnodes+i]); - Ke->values[(2*i+1)*numdof+2*j+1] += gauss->weight*Jdet*viscosity*(vdbasis[0*vnumnodes+j]*vdbasis[0*vnumnodes+i] + 2.*vdbasis[1*vnumnodes+j]*vdbasis[1*vnumnodes+i]); - } - for(int k=0;kvalues[(2*i+0)*numdof+2*vnumnodes+k] += gauss->weight*Jdet*FSreconditioning*(-pbasis[k]*vdbasis[0*vnumnodes+i]); - Ke->values[(2*i+1)*numdof+2*vnumnodes+k] += gauss->weight*Jdet*FSreconditioning*(-pbasis[k]*vdbasis[1*vnumnodes+i]); - } - } - /*Incompressibility*/ - for(int k=0;kvalues[(2*vnumnodes+k)*numdof+2*j+0] += gauss->weight*Jdet*(-FSreconditioning*vdbasis[0*vnumnodes+j]*pbasis[k]); - Ke->values[(2*vnumnodes+k)*numdof+2*j+1] += gauss->weight*Jdet*(-FSreconditioning*vdbasis[1*vnumnodes+j]*pbasis[k]); + for (int p=0;pvalues[(dim*vnumnodes+k)*numdof+dim*j+p] += gauss->weight*Jdet*(-FSreconditioning*vdbasis[p*vnumnodes+j]*pbasis[k]); + } } } @@ -3730,4 +3720,277 @@ xDelete(pbasis); return Ke; +}/*}}}*/ + +ElementMatrix* StressbalanceAnalysis::CreateKMatrixFSViscousGLS(Element* element){/*{{{*/ + + /*Intermediaries*/ + int i,dim; + IssmDouble viscosity,FSreconditioning,Jdet; + IssmDouble *xyz_list = NULL; + + /*Get problem dimension*/ + element->FindParam(&dim,DomainDimensionEnum); + + /*Fetch number of nodes and dof for this finite element*/ + int vnumnodes = element->NumberofNodesVelocity(); + int pnumnodes = element->NumberofNodesPressure(); + int numdof = vnumnodes*dim + pnumnodes; + + /* only work for P1 elements */ + if (dim == 2) {_assert_(vnumnodes==3);} + else if (dim == 3) {_assert_(vnumnodes==6);} + else {_error_("GLS is not implemented except for 2D and 3D problems.");} + + /*Prepare coordinate system list*/ + int* cs_list = xNew(vnumnodes+pnumnodes); + if(dim==2) for(i=0;iNewElementMatrix(FSvelocityEnum); + IssmDouble* vdbasis = xNew(dim*vnumnodes); + IssmDouble* pbasis = xNew(pnumnodes); + + /*Retrieve all inputs and parameters*/ + element->GetVerticesCoordinates(&xyz_list); + element->FindParam(&FSreconditioning,StressbalanceFSreconditioningEnum); + Input* vx_input=element->GetInput(VxEnum); _assert_(vx_input); + Input* vy_input=element->GetInput(VyEnum); _assert_(vy_input); + Input* vz_input = NULL; + if(dim==3){vz_input=element->GetInput(VzEnum); _assert_(vz_input);} + + /* prepare viscosity gradient for GLS */ + IssmDouble NodalViscosity[6]; + IssmDouble gradViscos[3]; + IssmDouble etapq, s, Tau, mk, hk; + IssmDouble hx, hy, hz; + IssmDouble SU[3*(3+1)*6]; + Gauss* vert = element->NewGauss(); + /* Compute the nodal values of the viscosity */ + for(int i=0;iGaussNode(element->element_type, i); + element->material->ViscosityFS(&NodalViscosity[i],dim,xyz_list,vert,vx_input,vy_input,vz_input); + } + delete vert; + + /* Start looping on the number of gaussian points: */ + Gauss* gauss = element->NewGauss(5); + while(gauss->next()){ + element->JacobianDeterminant(&Jdet,xyz_list,gauss); + element->NodalFunctionsDerivativesVelocity(vdbasis,xyz_list,gauss); + element->NodalFunctionsPressure(pbasis,gauss); + element->material->ViscosityFS(&viscosity,dim,xyz_list,gauss,vx_input,vy_input,vz_input); + + /* compute viscosity gradient at the gaussian point */ + element->ValueP1DerivativesOnGauss(&gradViscos[0],NodalViscosity,xyz_list,gauss); + + /* weight*detJ */ + s = gauss->weight*Jdet; + + /* GLS Stabilization */ + mk = 1.0 /3.0; + element->ElementSizes(&hx,&hy,&hz); + hk = max(max(hx, hy), hz); + // if(!element->IsOnDirichletBoundary()) { + Tau = - mk * hk * hk * 0.125 / viscosity; + + + for (int q=0;qvalues[(dim*i+p)*numdof+dim*j+p] += s*viscosity*(vdbasis[q*vnumnodes+j]*vdbasis[q*vnumnodes+i]); + /* All the entries */ + Ke->values[(dim*i+p)*numdof+dim*j+q] += s*viscosity*(vdbasis[p*vnumnodes+j]*vdbasis[q*vnumnodes+i]); + } + } + } + for(int k=0;kvalues[(dim*i+p)*numdof+dim*vnumnodes+k] += s*FSreconditioning*(-pbasis[k]*vdbasis[p*vnumnodes+i]); + } + } + } + /*Incompressibility*/ + for(int k=0;kvalues[(dim*vnumnodes+k)*numdof+dim*j+p] += s*(-FSreconditioning*vdbasis[p*vnumnodes+j]*pbasis[k]); + } + } + } + + /* GLS */ + for(int i=0;ivalues[i*numdof+j] += Tau * s * SU[p*numdof+i]*SU[p*numdof+j]; + } + } + } + } + } + + /*Transform Coordinate System*/ + element->TransformStiffnessMatrixCoord(Ke,cs_list); + + /*Clean up and return*/ + delete gauss; + xDelete(xyz_list); + xDelete(pbasis); + xDelete(vdbasis); + xDelete(cs_list); + return Ke; +}/*}}}*/ +ElementVector* StressbalanceAnalysis::CreatePVectorFSViscousGLS(Element* element){/*{{{*/ + + int i,dim; + IssmDouble Jdet,forcex,forcey,forcez; + IssmDouble *xyz_list = NULL; + + /*Get problem dimension*/ + element->FindParam(&dim,DomainDimensionEnum); + + /*Fetch number of nodes and dof for this finite element*/ + int vnumnodes = element->NumberofNodesVelocity(); + int pnumnodes = element->NumberofNodesPressure(); + + /*Prepare coordinate system list*/ + int* cs_list = xNew(vnumnodes+pnumnodes); + if(dim==2) for(i=0;iNewElementVector(FSvelocityEnum); + IssmDouble* vbasis = xNew(vnumnodes); + IssmDouble* vdbasis = xNew(dim*vnumnodes); + + /*Retrieve all inputs and parameters*/ + element->GetVerticesCoordinates(&xyz_list); + IssmDouble rho_ice =element->FindParam(MaterialsRhoIceEnum); + IssmDouble gravity =element->FindParam(ConstantsGEnum); + Input* loadingforcex_input=element->GetInput(LoadingforceXEnum); _assert_(loadingforcex_input); + Input* loadingforcey_input=element->GetInput(LoadingforceYEnum); _assert_(loadingforcey_input); + Input* loadingforcez_input=NULL; + if(dim==3){ + loadingforcez_input=element->GetInput(LoadingforceZEnum); _assert_(loadingforcez_input); + } + + /* prepare viscosity gradient for GLS */ + Input* vx_input=element->GetInput(VxEnum); _assert_(vx_input); + Input* vy_input=element->GetInput(VyEnum); _assert_(vy_input); + Input* vz_input = NULL; + if(dim==3){vz_input=element->GetInput(VzEnum); _assert_(vz_input);} + + IssmDouble viscosity,FSreconditioning; + IssmDouble NodalViscosity[6]; + IssmDouble gradViscos[3]; + IssmDouble etapq, s, Tau, mk, hk; + IssmDouble hx, hy, hz; + IssmDouble SU[3*(3+1)*6]; + Gauss* vert = element->NewGauss(); + + element->FindParam(&FSreconditioning,StressbalanceFSreconditioningEnum); + /* Compute the nodal values of the viscosity */ + for(int i=0;iGaussNode(element->element_type, i); + element->material->ViscosityFS(&NodalViscosity[i],dim,xyz_list,vert,vx_input,vy_input,vz_input); + } + delete vert; + + /* Start looping on the number of gaussian points: */ + Gauss* gauss=element->NewGauss(5); + while(gauss->next()){ + + element->JacobianDeterminant(&Jdet,xyz_list,gauss); + element->NodalFunctionsVelocity(vbasis,gauss); + element->NodalFunctionsDerivativesVelocity(vdbasis,xyz_list,gauss); + + loadingforcex_input->GetInputValue(&forcex,gauss); + loadingforcey_input->GetInputValue(&forcey,gauss); + if(dim==3) loadingforcez_input->GetInputValue(&forcez,gauss); + + /* compute viscosity gradient at the gaussian point */ + element->ValueP1DerivativesOnGauss(&gradViscos[0],NodalViscosity,xyz_list,gauss); + /* weight*detJ */ + s = gauss->weight*Jdet; + /* GLS Stabilization */ + element->material->ViscosityFS(&viscosity,dim,xyz_list,gauss,vx_input,vy_input,vz_input); + mk = 1.0 /3.0; + element->ElementSizes(&hx,&hy,&hz); + hk = max(max(hx, hy), hz); + // if(!element->IsOnDirichletBoundary()) { + Tau = - mk * hk * hk * 0.125 / viscosity; + // } + // else { + // Tau = 0.0; + // } + + for (int q=0;qvalues[i*dim+0] += +rho_ice*forcex *Jdet*gauss->weight*vbasis[i]; + pe->values[i*dim+1] += +rho_ice*forcey *Jdet*gauss->weight*vbasis[i]; + if(dim==3) pe->values[i*dim+2] += +rho_ice*forcez*Jdet*gauss->weight*vbasis[i]; + + pe->values[i*dim+dim-1] += -rho_ice*gravity*Jdet*gauss->weight*vbasis[i]; + } + + for(int i=0;ivalues[i] += Tau*Jdet*gauss->weight*rho_ice*(-gravity)*SU[(dim-1)*(dim+1)*vnumnodes+i]; + } + } + + /*Transform coordinate system*/ + element->TransformLoadVectorCoord(pe,cs_list); + + /*Clean up and return*/ + delete gauss; + xDelete(cs_list); + xDelete(vbasis); + xDelete(vdbasis); + xDelete(xyz_list); + return pe; }/*}}}*/ @@ -4335,4 +4598,13 @@ delete pe4; } + else if(fe_FS==P1P1GLSEnum) { + ElementVector* pe1=CreatePVectorFSViscousGLS(element); + ElementVector* pe2=CreatePVectorFSShelf(element); + ElementVector* pe3=CreatePVectorFSFront(element); + pe =new ElementVector(pe1,pe2,pe3); + delete pe1; + delete pe2; + delete pe3; + } else{ ElementVector* pe1=CreatePVectorFSViscous(element); @@ -4396,11 +4668,7 @@ pe->values[i*dim+0] += +rho_ice*forcex *Jdet*gauss->weight*vbasis[i]; pe->values[i*dim+1] += +rho_ice*forcey *Jdet*gauss->weight*vbasis[i]; - if(dim==3){ - pe->values[i*dim+2] += +rho_ice*forcez*Jdet*gauss->weight*vbasis[i]; - pe->values[i*dim+2] += -rho_ice*gravity*Jdet*gauss->weight*vbasis[i]; - } - else{ - pe->values[i*dim+1] += -rho_ice*gravity*Jdet*gauss->weight*vbasis[i]; - } + if(dim==3) pe->values[i*dim+2] += +rho_ice*forcez*Jdet*gauss->weight*vbasis[i]; + + pe->values[i*dim+dim-1] += -rho_ice*gravity*Jdet*gauss->weight*vbasis[i]; } } @@ -4415,4 +4683,190 @@ xDelete(xyz_list); return pe; +}/*}}}*/ +ElementMatrix* StressbalanceAnalysis::CreateKMatrixFSFrictionNitsche(Element* element){/*{{{*/ + + if(element->IsFloating() || !element->IsOnBase()) return NULL; + + /*If on water or not FS, skip stiffness: */ + int approximation; + element->GetInputValue(&approximation,ApproximationEnum); + if(approximation!=FSApproximationEnum && approximation!=SSAFSApproximationEnum && approximation!=HOFSApproximationEnum) return NULL; + + /*Intermediaries*/ + bool mainlyfloating; + int dim,domaintype; + int friction_style,point1; + IssmDouble alpha2,Jdet,fraction1,fraction2; + IssmDouble viscosity, FSreconditioning; + IssmDouble gllevelset,phi=1.; + IssmDouble *xyz_list_base = NULL; + IssmDouble *xyz_list = NULL; + Gauss* gauss = NULL; + /*coefficient of Nitsche's method*/ + IssmDouble gamma, hx, hy, hz; + IssmDouble pnu, qnv; + + /*Get problem dimension*/ + element->FindParam(&dim,DomainDimensionEnum); + element->FindParam(&domaintype,DomainTypeEnum); + element->FindParam(&FSreconditioning,StressbalanceFSreconditioningEnum); + + /*Fetch number of nodes and dof for this finite element*/ + int vnumnodes = element->NumberofNodesVelocity(); + int pnumnodes = element->NumberofNodesPressure(); + int numdof = vnumnodes*dim + pnumnodes; + + /*Initialize Element matrix and vectors*/ + ElementMatrix* Ke = element->NewElementMatrix(FSvelocityEnum); + IssmDouble* vbasis = xNew(vnumnodes); + IssmDouble* vdbasis = xNew(dim*vnumnodes); + IssmDouble* pbasis = xNew(pnumnodes); + IssmDouble* vdbasisn = xNew(vnumnodes); + IssmDouble* nsigma = xNew(numdof); + IssmDouble bed_normal[3], bed_tangent[6]; + + /*Prepare coordinate system list*/ + int* cs_list = xNew(vnumnodes+pnumnodes); + if(dim==2) for(int i=0;iGetVerticesCoordinatesBase(&xyz_list_base); + element->GetVerticesCoordinates(&xyz_list); + element->FindParam(&friction_style,GroundinglineFrictionInterpolationEnum); + Input* gllevelset_input = NULL; + Input* vx_input=element->GetInput(VxEnum); _assert_(vx_input); + Input* vy_input=element->GetInput(VyEnum); _assert_(vy_input); + Input* vz_input = NULL; + if(dim==3){vz_input=element->GetInput(VzEnum); _assert_(vz_input);} + + /*build friction object, used later on: */ + Friction* friction=new Friction(element,dim==3?3:1); + + /*Recover portion of element that is grounded*/ + if(!(friction_style==SubelementFriction2Enum)) phi=element->GetGroundedPortion(xyz_list_base); + if(friction_style==SubelementFriction2Enum){ + if(domaintype==Domain2DverticalEnum) _error_("Subelement Friction 2 not implemented yet for Flowline"); + gllevelset_input=element->GetInput(MaskOceanLevelsetEnum); _assert_(gllevelset_input); + element->GetGroundedPart(&point1,&fraction1,&fraction2,&mainlyfloating); + gauss=element->NewGaussBase(3); + } + else{ + gauss=element->NewGaussBase(3); + } + + element->NormalBase(&bed_normal[0],xyz_list_base); + element->TangentBase(&bed_tangent[0], &bed_normal[0]); + element->ElementSizes(&hx,&hy,&hz); + element->FindParam(&gamma,FeFSNitscheGammaEnum); + + gamma = gamma / sqrt(hx*hx+hy*hy+hz*hz); + + /* Start looping on the number of gaussian points: */ + while(gauss->next()){ + friction->GetAlpha2(&alpha2,gauss); + if(friction_style==SubelementFriction1Enum) alpha2=phi*alpha2; + else if(friction_style==SubelementFriction2Enum){ + gllevelset_input->GetInputValue(&gllevelset, gauss); + if(gllevelset<0.) alpha2=0.; + } + else if(friction_style==NoFrictionOnPartiallyFloatingEnum){ + if (phi<0.99999999) alpha2=0.; + } + else _error_("friction interpolation "<JacobianDeterminantBase(&Jdet,xyz_list_base,gauss); + element->NodalFunctionsVelocity(vbasis,gauss); + /* The full element is needed for calculating derivative of the basis functions */ + element->NodalFunctionsDerivativesVelocity(vdbasis,xyz_list,gauss); + element->NodalFunctionsPressure(pbasis,gauss); + element->material->ViscosityFS(&viscosity,dim,xyz_list,gauss,vx_input,vy_input,vz_input); + + /*Frictions*/ + for(int i=0;ivalues[(dim*i+p)*numdof+dim*j+q] += alpha2*gauss->weight*Jdet*vbasis[i]*vbasis[j]*bed_tangent[d*3+p]*bed_tangent[d*3+q]; + } + } + } + } + } + /* -------- Nitsche terms -------- */ + /* n\cdot\nabla\phi*/ + for(int i=0;ivalues[(dim*i+p)*numdof+dim*j+q] += gauss->weight*Jdet * (-gamma) * bed_normal[q] * vbasis[j]*bed_normal[p] * vbasis[i]; + /* sigma(u)*(n*v) */ + Ke->values[(dim*i+p)*numdof+dim*j+q] += (- gauss->weight*Jdet)* nsigma[p*vnumnodes+i] * bed_normal[q] * vbasis[j]; + /* sigma(v)*(n*u) */ + Ke->values[(dim*i+p)*numdof+dim*j+q] += (- gauss->weight*Jdet)* bed_normal[p] * vbasis[i] * nsigma[q*vnumnodes+j]; + } + } + } + } + + /* pressure x velocity component A12 */ + for(int k=0;kvalues[(dim*i+p)*numdof+dim*vnumnodes+k] += gauss->weight*Jdet * pnu * (- bed_normal[p] * vbasis[i]); + } + } + } + /* velocity x pressure component A21 */ + for(int k=0;kvalues[(dim*vnumnodes+k)*numdof+dim*j+q] += gauss->weight*Jdet * ( - bed_normal[q] * vbasis[j]) * qnv; + } + } + } + } + + /*Transform Coordinate System*/ + element->TransformStiffnessMatrixCoord(Ke,cs_list); + + /*Clean up and return*/ + delete gauss; + delete friction; + xDelete(xyz_list_base); + xDelete(xyz_list); + xDelete(vbasis); + xDelete(vdbasis); + xDelete(pbasis); + xDelete(vdbasisn); + xDelete(nsigma); + xDelete(cs_list); + return Ke; }/*}}}*/ #endif Index: /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.h =================================================================== --- /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.h (revision 25513) +++ /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.h (revision 25514) @@ -73,6 +73,8 @@ ElementMatrix* CreateKMatrixFS(Element* element); ElementMatrix* CreateKMatrixFSFriction(Element* element); + ElementMatrix* CreateKMatrixFSFrictionNitsche(Element* element); ElementMatrix* CreateKMatrixFSShelf(Element* element); ElementMatrix* CreateKMatrixFSViscous(Element* element); + ElementMatrix* CreateKMatrixFSViscousGLS(Element* element); ElementMatrix* CreateKMatrixFSViscousLA(Element* element); ElementMatrix* CreateKMatrixFSViscousXTH(Element* element); @@ -85,4 +87,5 @@ ElementVector* CreatePVectorFSStress(Element* element); ElementVector* CreatePVectorFSViscous(Element* element); + ElementVector* CreatePVectorFSViscousGLS(Element* element); ElementVector* CreatePVectorFSViscousLA(Element* element); ElementVector* CreatePVectorFSViscousXTH(Element* element); Index: /issm/trunk-jpl/src/c/classes/Elements/Element.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Element.h (revision 25513) +++ /issm/trunk-jpl/src/c/classes/Elements/Element.h (revision 25514) @@ -343,4 +343,5 @@ virtual void StressIntensityFactor(void)=0; virtual IssmDouble SurfaceArea(void)=0; + virtual void TangentBase(IssmDouble* bed_tangent,IssmDouble* bed_normal){_error_("not implemented yet");}; virtual int TensorInterpolation()=0; virtual IssmDouble TimeAdapt()=0; Index: /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp (revision 25513) +++ /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp (revision 25514) @@ -3679,4 +3679,44 @@ return S; } +} +/*}}}*/ +void Penta::TangentBase(IssmDouble* bed_tangent,IssmDouble* bed_normal){/*{{{*/ + /* + To compute the two tangent bed_tangent[0:2] and bed_tangent[3:5] from the given normal vecotr. + */ + IssmDouble n1, n2, n3; + IssmDouble tangent_norm; + + n1 = fabs(bed_normal[0]); + n2 = fabs(bed_normal[1]); + n3 = fabs(bed_normal[2]); + + /* For the first tangent, on x-y plane in most cases*/ + if ((n1<=n3) && (n2<=n3)) { + bed_tangent[0] = 0.0; + bed_tangent[1] = -bed_normal[2]; + bed_tangent[2] = bed_normal[1]; + } + else { + bed_tangent[0] = -bed_normal[1]; + bed_tangent[1] = bed_normal[0]; + bed_tangent[2] = 0.0; + } + tangent_norm = sqrt(bed_tangent[0]*bed_tangent[0]+bed_tangent[1]*bed_tangent[1]+bed_tangent[2]*bed_tangent[2]); + for(int i=0;i<3;i++) bed_tangent[i] = bed_tangent[i]/tangent_norm; + /* The second tangent*/ + bed_tangent[3] = bed_normal[1]*bed_tangent[2] - bed_normal[2]*bed_tangent[1]; + bed_tangent[4] = bed_normal[2]*bed_tangent[0] - bed_normal[0]*bed_tangent[2]; + bed_tangent[5] = bed_normal[0]*bed_tangent[1] - bed_normal[1]*bed_tangent[0]; + tangent_norm = sqrt(bed_tangent[3]*bed_tangent[3]+bed_tangent[4]*bed_tangent[4]+bed_tangent[5]*bed_tangent[5]); + for(int i=3;i<6;i++) bed_tangent[i] = bed_tangent[i]/tangent_norm; + + IssmDouble checksum = 0.0; + for (int i=0;i<3;i++) { + checksum += bed_normal[i]*bed_tangent[i]; + checksum += bed_normal[i]*bed_tangent[i+3]; + checksum += bed_tangent[i]*bed_tangent[i+3]; + } + _assert_(fabs(checksum)<1e-10); } /*}}}*/ Index: /issm/trunk-jpl/src/c/classes/Elements/Penta.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Penta.h (revision 25513) +++ /issm/trunk-jpl/src/c/classes/Elements/Penta.h (revision 25514) @@ -178,4 +178,5 @@ void StrainRateperpendicular(); IssmDouble SurfaceArea(void); + void TangentBase(IssmDouble* bed_tangent,IssmDouble* bed_normal); int TensorInterpolation(){_error_("not implemented yet");}; IssmDouble TimeAdapt(); Index: /issm/trunk-jpl/src/c/classes/gauss/GaussPenta.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/gauss/GaussPenta.cpp (revision 25513) +++ /issm/trunk-jpl/src/c/classes/gauss/GaussPenta.cpp (revision 25514) @@ -571,4 +571,5 @@ switch(finiteelement){ case P1Enum: case P1DGEnum: + case P1P1GLSEnum: case P1P1Enum: /* added to allow P1-P1 GLS */ switch(iv){ case 0: coord1=1.; coord2=0.; coord3=0.; coord4=-1.; break; Index: /issm/trunk-jpl/src/c/classes/gauss/GaussTria.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/gauss/GaussTria.cpp (revision 25513) +++ /issm/trunk-jpl/src/c/classes/gauss/GaussTria.cpp (revision 25514) @@ -567,4 +567,5 @@ break; case P1Enum: case P1DGEnum: + case P1P1GLSEnum: case P1P1Enum:/* added to allow P1-P1 GLS */ switch(iv){ case 0: coord1=1.; coord2=0.; coord3=0.; break; Index: /issm/trunk-jpl/src/c/cores/stressbalance_core.cpp =================================================================== --- /issm/trunk-jpl/src/c/cores/stressbalance_core.cpp (revision 25513) +++ /issm/trunk-jpl/src/c/cores/stressbalance_core.cpp (revision 25514) @@ -19,5 +19,5 @@ bool dakota_analysis,control_analysis; int domaintype; - bool isSIA,isSSA,isL1L2,isHO,isFS; + bool isSIA,isSSA,isL1L2,isHO,isFS,isNitsche; bool save_results; int solution_type; @@ -33,4 +33,5 @@ femmodel->parameters->FindParam(&isHO,FlowequationIsHOEnum); femmodel->parameters->FindParam(&isFS,FlowequationIsFSEnum); + femmodel->parameters->FindParam(&isNitsche,FlowequationIsNitscheEnum); femmodel->parameters->FindParam(&save_results,SaveResultsEnum); femmodel->parameters->FindParam(&solution_type,SolutionTypeEnum); @@ -46,5 +47,7 @@ bedslope_core(femmodel); femmodel->SetCurrentConfiguration(StressbalanceAnalysisEnum); - ResetFSBasalBoundaryConditionx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters); + if (!isNitsche) { + ResetFSBasalBoundaryConditionx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters); + } /*We need basal melt rates for the shelf dampening*/ Index: /issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateNodes.cpp =================================================================== --- /issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateNodes.cpp (revision 25513) +++ /issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateNodes.cpp (revision 25514) @@ -66,4 +66,8 @@ numnodes = iomodel->numberofvertices+iomodel->numberofedges; break; + case P1P1Enum: case P1P1GLSEnum: + /*P1 velocity + P1 pressure element with GLS stabilization*/ + numnodes = (iomodel->numberofvertices) + iomodel->numberofvertices; + break; case MINIEnum: case MINIcondensedEnum: /*P1+ velocity (bubble statically condensed), P1 pressure*/ @@ -127,4 +131,8 @@ FacesPartitioning(iomodel); numnodes = iomodel->numberofvertices+iomodel->numberofedges+2*iomodel->numberofverticaledges+3*iomodel->numberofverticalfaces; + break; + case P1P1Enum: case P1P1GLSEnum: + /*P1 velocity + P1 pressure element with GLS stabilization*/ + numnodes = (iomodel->numberofvertices) + iomodel->numberofvertices; break; case MINIEnum: case MINIcondensedEnum: @@ -207,4 +215,15 @@ element_node_ids[4]=iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[3*i+1]; element_node_ids[5]=iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[3*i+2]; + break; + case P1P1Enum: case P1P1GLSEnum: + element_numnodes = 3+3; + element_node_ids[0]=iomodel->elements[3*i+0]-1; + element_node_ids[1]=iomodel->elements[3*i+1]-1; + element_node_ids[2]=iomodel->elements[3*i+2]-1; + for(int n=0;n<3;n++) nodes_approx[element_node_ids[n]] = FSvelocityEnum; + element_node_ids[3]=iomodel->numberofvertices+iomodel->elements[3*i+0]-1; + element_node_ids[4]=iomodel->numberofvertices+iomodel->elements[3*i+1]-1; + element_node_ids[5]=iomodel->numberofvertices+iomodel->elements[3*i+2]-1; + for(int n=3;n<6;n++) nodes_approx[element_node_ids[n]] = FSpressureEnum; break; case MINIEnum: case MINIcondensedEnum: @@ -408,4 +427,23 @@ element_node_ids[28]=iomodel->numberofvertices+iomodel->numberofedges+2*iomodel->numberofverticaledges+3*iomodel->elementtoverticalfaceconnectivity[3*i+1]+2; /* 3/4 vertical face 2*/ element_node_ids[29]=iomodel->numberofvertices+iomodel->numberofedges+2*iomodel->numberofverticaledges+3*iomodel->elementtoverticalfaceconnectivity[3*i+2]+2; /* 3/4 vertical face 3*/ + break; + case P1P1Enum: case P1P1GLSEnum: + /*P1-P1 very similar to MINI, but no DOF on the element + TODO: add if(!approximations) or not?*/ + element_numnodes = 6+6; + element_node_ids[ 0]=iomodel->elements[6*i+0]-1; + element_node_ids[ 1]=iomodel->elements[6*i+1]-1; + element_node_ids[ 2]=iomodel->elements[6*i+2]-1; + element_node_ids[ 3]=iomodel->elements[6*i+3]-1; + element_node_ids[ 4]=iomodel->elements[6*i+4]-1; + element_node_ids[ 5]=iomodel->elements[6*i+5]-1; + if(!approximations) for(int n=0;n<6;n ++) nodes_approx[element_node_ids[n]] = FSvelocityEnum; + element_node_ids[ 6]=iomodel->numberofvertices+iomodel->elements[6*i+0]-1; + element_node_ids[ 7]=iomodel->numberofvertices+iomodel->elements[6*i+1]-1; + element_node_ids[ 8]=iomodel->numberofvertices+iomodel->elements[6*i+2]-1; + element_node_ids[ 9]=iomodel->numberofvertices+iomodel->elements[6*i+3]-1; + element_node_ids[10]=iomodel->numberofvertices+iomodel->elements[6*i+4]-1; + element_node_ids[11]=iomodel->numberofvertices+iomodel->elements[6*i+5]-1; + if(!approximations) for(int n=6;n<12;n++) nodes_approx[element_node_ids[n]] = FSpressureEnum; break; case MINIEnum: case MINIcondensedEnum: Index: /issm/trunk-jpl/src/c/modules/SetActiveNodesLSMx/SetActiveNodesLSMx.cpp =================================================================== --- /issm/trunk-jpl/src/c/modules/SetActiveNodesLSMx/SetActiveNodesLSMx.cpp (revision 25513) +++ /issm/trunk-jpl/src/c/modules/SetActiveNodesLSMx/SetActiveNodesLSMx.cpp (revision 25514) @@ -24,4 +24,5 @@ /*include switch for elements with multiple different sets of nodes*/ switch(element->GetElementType()){ + case P1P1GLSEnum: case P1P1Enum:/* added to allow P1-P1 GLS */ case MINIEnum:case MINIcondensedEnum: case TaylorHoodEnum:case XTaylorHoodEnum:case LATaylorHoodEnum: Index: /issm/trunk-jpl/src/c/shared/Enum/Enum.vim =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/Enum.vim (revision 25513) +++ /issm/trunk-jpl/src/c/shared/Enum/Enum.vim (revision 25514) @@ -148,4 +148,6 @@ syn keyword cConstant FlowequationIsSIAEnum syn keyword cConstant FlowequationIsSSAEnum +syn keyword cConstant FlowequationIsNitscheEnum +syn keyword cConstant FeFSNitscheGammaEnum syn keyword cConstant FrictionCouplingEnum syn keyword cConstant FrictionDeltaEnum Index: /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h (revision 25513) +++ /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h (revision 25514) @@ -142,4 +142,6 @@ FlowequationIsSIAEnum, FlowequationIsSSAEnum, + FlowequationIsNitscheEnum, + FeFSNitscheGammaEnum, FrictionCouplingEnum, FrictionDeltaEnum, Index: /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp (revision 25513) +++ /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp (revision 25514) @@ -150,4 +150,6 @@ case FlowequationIsSIAEnum : return "FlowequationIsSIA"; case FlowequationIsSSAEnum : return "FlowequationIsSSA"; + case FlowequationIsNitscheEnum : return "FlowequationIsNitsche"; + case FeFSNitscheGammaEnum : return "FeFSNitscheGamma"; case FrictionCouplingEnum : return "FrictionCoupling"; case FrictionDeltaEnum : return "FrictionDelta"; Index: /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp (revision 25513) +++ /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp (revision 25514) @@ -153,4 +153,6 @@ else if (strcmp(name,"FlowequationIsSIA")==0) return FlowequationIsSIAEnum; else if (strcmp(name,"FlowequationIsSSA")==0) return FlowequationIsSSAEnum; + else if (strcmp(name,"FlowequationIsNitsche")==0) return FlowequationIsNitscheEnum; + else if (strcmp(name,"FeFSNitscheGamma")==0) return FeFSNitscheGammaEnum; else if (strcmp(name,"FrictionCoupling")==0) return FrictionCouplingEnum; else if (strcmp(name,"FrictionDelta")==0) return FrictionDeltaEnum; @@ -258,10 +260,10 @@ else if (strcmp(name,"MasstransportNumRequestedOutputs")==0) return MasstransportNumRequestedOutputsEnum; else if (strcmp(name,"MasstransportPenaltyFactor")==0) return MasstransportPenaltyFactorEnum; - else if (strcmp(name,"MasstransportRequestedOutputs")==0) return MasstransportRequestedOutputsEnum; - else if (strcmp(name,"MasstransportStabilization")==0) return MasstransportStabilizationEnum; else stage=3; } if(stage==3){ - if (strcmp(name,"MaterialsBeta")==0) return MaterialsBetaEnum; + if (strcmp(name,"MasstransportRequestedOutputs")==0) return MasstransportRequestedOutputsEnum; + else if (strcmp(name,"MasstransportStabilization")==0) return MasstransportStabilizationEnum; + else if (strcmp(name,"MaterialsBeta")==0) return MaterialsBetaEnum; else if (strcmp(name,"MaterialsEarthDensity")==0) return MaterialsEarthDensityEnum; else if (strcmp(name,"MaterialsEffectiveconductivityAveraging")==0) return MaterialsEffectiveconductivityAveragingEnum; @@ -381,10 +383,10 @@ else if (strcmp(name,"SmbIsdelta18o")==0) return SmbIsdelta18oEnum; else if (strcmp(name,"SmbIsdensification")==0) return SmbIsdensificationEnum; - else if (strcmp(name,"SmbIsfirnwarming")==0) return SmbIsfirnwarmingEnum; - else if (strcmp(name,"SmbIsgraingrowth")==0) return SmbIsgraingrowthEnum; else stage=4; } if(stage==4){ - if (strcmp(name,"SmbIsmelt")==0) return SmbIsmeltEnum; + if (strcmp(name,"SmbIsfirnwarming")==0) return SmbIsfirnwarmingEnum; + else if (strcmp(name,"SmbIsgraingrowth")==0) return SmbIsgraingrowthEnum; + else if (strcmp(name,"SmbIsmelt")==0) return SmbIsmeltEnum; else if (strcmp(name,"SmbIsmungsm")==0) return SmbIsmungsmEnum; else if (strcmp(name,"SmbIsprecipscaled")==0) return SmbIsprecipscaledEnum; @@ -504,10 +506,10 @@ else if (strcmp(name,"BasalforcingsPicoBoxId")==0) return BasalforcingsPicoBoxIdEnum; else if (strcmp(name,"BasalforcingsPicoOverturningCoeff")==0) return BasalforcingsPicoOverturningCoeffEnum; - else if (strcmp(name,"BasalforcingsPicoSubShelfOceanOverturning")==0) return BasalforcingsPicoSubShelfOceanOverturningEnum; - else if (strcmp(name,"BasalforcingsPicoSubShelfOceanSalinity")==0) return BasalforcingsPicoSubShelfOceanSalinityEnum; else stage=5; } if(stage==5){ - if (strcmp(name,"BasalforcingsPicoSubShelfOceanTemp")==0) return BasalforcingsPicoSubShelfOceanTempEnum; + if (strcmp(name,"BasalforcingsPicoSubShelfOceanOverturning")==0) return BasalforcingsPicoSubShelfOceanOverturningEnum; + else if (strcmp(name,"BasalforcingsPicoSubShelfOceanSalinity")==0) return BasalforcingsPicoSubShelfOceanSalinityEnum; + else if (strcmp(name,"BasalforcingsPicoSubShelfOceanTemp")==0) return BasalforcingsPicoSubShelfOceanTempEnum; else if (strcmp(name,"BasalStressx")==0) return BasalStressxEnum; else if (strcmp(name,"BasalStressy")==0) return BasalStressyEnum; @@ -627,10 +629,10 @@ else if (strcmp(name,"HydrologydcMaskEplactiveElt")==0) return HydrologydcMaskEplactiveEltEnum; else if (strcmp(name,"HydrologydcMaskEplactiveNode")==0) return HydrologydcMaskEplactiveNodeEnum; - else if (strcmp(name,"HydrologydcMaskThawedElt")==0) return HydrologydcMaskThawedEltEnum; - else if (strcmp(name,"HydrologydcMaskThawedNode")==0) return HydrologydcMaskThawedNodeEnum; else stage=6; } if(stage==6){ - if (strcmp(name,"HydrologydcSedimentTransmitivity")==0) return HydrologydcSedimentTransmitivityEnum; + if (strcmp(name,"HydrologydcMaskThawedElt")==0) return HydrologydcMaskThawedEltEnum; + else if (strcmp(name,"HydrologydcMaskThawedNode")==0) return HydrologydcMaskThawedNodeEnum; + else if (strcmp(name,"HydrologydcSedimentTransmitivity")==0) return HydrologydcSedimentTransmitivityEnum; else if (strcmp(name,"HydrologyDrainageRate")==0) return HydrologyDrainageRateEnum; else if (strcmp(name,"HydrologyEnglacialInput")==0) return HydrologyEnglacialInputEnum; @@ -750,10 +752,10 @@ else if (strcmp(name,"SmbDswrf")==0) return SmbDswrfEnum; else if (strcmp(name,"SmbDzAdd")==0) return SmbDzAddEnum; - else if (strcmp(name,"SmbDz")==0) return SmbDzEnum; - else if (strcmp(name,"SmbDzMin")==0) return SmbDzMinEnum; else stage=7; } if(stage==7){ - if (strcmp(name,"SmbDzTop")==0) return SmbDzTopEnum; + if (strcmp(name,"SmbDz")==0) return SmbDzEnum; + else if (strcmp(name,"SmbDzMin")==0) return SmbDzMinEnum; + else if (strcmp(name,"SmbDzTop")==0) return SmbDzTopEnum; else if (strcmp(name,"SmbDzini")==0) return SmbDziniEnum; else if (strcmp(name,"SmbEAir")==0) return SmbEAirEnum; @@ -873,10 +875,10 @@ else if (strcmp(name,"VxMesh")==0) return VxMeshEnum; else if (strcmp(name,"VxObs")==0) return VxObsEnum; - else if (strcmp(name,"VyAverage")==0) return VyAverageEnum; - else if (strcmp(name,"Vy")==0) return VyEnum; else stage=8; } if(stage==8){ - if (strcmp(name,"VyMesh")==0) return VyMeshEnum; + if (strcmp(name,"VyAverage")==0) return VyAverageEnum; + else if (strcmp(name,"Vy")==0) return VyEnum; + else if (strcmp(name,"VyMesh")==0) return VyMeshEnum; else if (strcmp(name,"VyObs")==0) return VyObsEnum; else if (strcmp(name,"Vz")==0) return VzEnum; @@ -996,10 +998,10 @@ else if (strcmp(name,"AdaptiveTimestepping")==0) return AdaptiveTimesteppingEnum; else if (strcmp(name,"AdjointBalancethickness2Analysis")==0) return AdjointBalancethickness2AnalysisEnum; - else if (strcmp(name,"AdjointBalancethicknessAnalysis")==0) return AdjointBalancethicknessAnalysisEnum; - else if (strcmp(name,"AdjointHorizAnalysis")==0) return AdjointHorizAnalysisEnum; else stage=9; } if(stage==9){ - if (strcmp(name,"AggressiveMigration")==0) return AggressiveMigrationEnum; + if (strcmp(name,"AdjointBalancethicknessAnalysis")==0) return AdjointBalancethicknessAnalysisEnum; + else if (strcmp(name,"AdjointHorizAnalysis")==0) return AdjointHorizAnalysisEnum; + else if (strcmp(name,"AggressiveMigration")==0) return AggressiveMigrationEnum; else if (strcmp(name,"AmrBamg")==0) return AmrBamgEnum; else if (strcmp(name,"AmrNeopz")==0) return AmrNeopzEnum; @@ -1119,10 +1121,10 @@ else if (strcmp(name,"Gset")==0) return GsetEnum; else if (strcmp(name,"Gsl")==0) return GslEnum; - else if (strcmp(name,"HOApproximation")==0) return HOApproximationEnum; - else if (strcmp(name,"HOFSApproximation")==0) return HOFSApproximationEnum; else stage=10; } if(stage==10){ - if (strcmp(name,"Hook")==0) return HookEnum; + if (strcmp(name,"HOApproximation")==0) return HOApproximationEnum; + else if (strcmp(name,"HOFSApproximation")==0) return HOFSApproximationEnum; + else if (strcmp(name,"Hook")==0) return HookEnum; else if (strcmp(name,"HydrologyDCEfficientAnalysis")==0) return HydrologyDCEfficientAnalysisEnum; else if (strcmp(name,"HydrologyDCInefficientAnalysis")==0) return HydrologyDCInefficientAnalysisEnum; @@ -1242,10 +1244,10 @@ else if (strcmp(name,"P1P1GLS")==0) return P1P1GLSEnum; else if (strcmp(name,"P1bubble")==0) return P1bubbleEnum; - else if (strcmp(name,"P1bubblecondensed")==0) return P1bubblecondensedEnum; - else if (strcmp(name,"P1xP2")==0) return P1xP2Enum; else stage=11; } if(stage==11){ - if (strcmp(name,"P1xP3")==0) return P1xP3Enum; + if (strcmp(name,"P1bubblecondensed")==0) return P1bubblecondensedEnum; + else if (strcmp(name,"P1xP2")==0) return P1xP2Enum; + else if (strcmp(name,"P1xP3")==0) return P1xP3Enum; else if (strcmp(name,"P1xP4")==0) return P1xP4Enum; else if (strcmp(name,"P2")==0) return P2Enum; @@ -1365,10 +1367,10 @@ else if (strcmp(name,"MeshZ")==0) return MeshZEnum; else if (strcmp(name,"NearestInterp")==0) return NearestInterpEnum; - else if (strcmp(name,"OutputdefinitionList")==0) return OutputdefinitionListEnum; - else if (strcmp(name,"SealevelObs")==0) return SealevelObsEnum; else stage=12; } if(stage==12){ - if (strcmp(name,"SealevelWeights")==0) return SealevelWeightsEnum; + if (strcmp(name,"OutputdefinitionList")==0) return OutputdefinitionListEnum; + else if (strcmp(name,"SealevelObs")==0) return SealevelObsEnum; + else if (strcmp(name,"SealevelWeights")==0) return SealevelWeightsEnum; else if (strcmp(name,"StrainRate")==0) return StrainRateEnum; else if (strcmp(name,"StressTensor")==0) return StressTensorEnum; Index: /issm/trunk-jpl/src/m/classes/flowequation.m =================================================================== --- /issm/trunk-jpl/src/m/classes/flowequation.m (revision 25513) +++ /issm/trunk-jpl/src/m/classes/flowequation.m (revision 25514) @@ -11,4 +11,6 @@ isHO = 0; isFS = 0; + isNitscheBC = 0; + FSNitscheGamma = 1e6; fe_SSA = ''; fe_HO = ''; @@ -98,4 +100,6 @@ md = checkfield(md,'fieldname','flowequation.isHO','numel',[1],'values',[0 1]); md = checkfield(md,'fieldname','flowequation.isFS','numel',[1],'values',[0 1]); + md = checkfield(md,'fieldname','flowequation.isNitscheBC','numel',[1],'values',[0 1]); + md = checkfield(md,'fieldname','flowequation.FSNitscheGamma','numel',[1], '>=', 0.); md = checkfield(md,'fieldname','flowequation.fe_SSA','values',{'P1','P1bubble','P1bubblecondensed','P2','P2bubble'}); md = checkfield(md,'fieldname','flowequation.fe_HO' ,'values',{'P1','P1bubble','P1bubblecondensed','P1xP2','P2xP1','P2','P2bubble','P1xP3','P1xP4','P2xP4'}); @@ -144,4 +148,6 @@ fielddisplay(self,'isHO','is the Higher-Order (HO) approximation used ?'); fielddisplay(self,'isFS','are the Full-FS (FS) equations used ?'); + fielddisplay(self,'isNitscheBC','is weakly imposed condition used?'); + fielddisplay(self,'FSNitscheGamma','Gamma value for the Nitsche term, by default gamma=1e6?'); fielddisplay(self,'fe_SSA','Finite Element for SSA ''P1'', ''P1bubble'' ''P1bubblecondensed'' ''P2'''); fielddisplay(self,'fe_HO', 'Finite Element for HO ''P1'' ''P1bubble'' ''P1bubblecondensed'' ''P1xP2'' ''P2xP1'' ''P2'''); @@ -160,4 +166,6 @@ WriteData(fid,prefix,'object',self,'fieldname','isHO','format','Boolean'); WriteData(fid,prefix,'object',self,'fieldname','isFS','format','Boolean'); + WriteData(fid,prefix,'object',self,'fieldname','isNitscheBC','format','Boolean'); + WriteData(fid,prefix,'object',self,'fieldname','FSNitscheGamma','format','Double'); WriteData(fid,prefix,'object',self,'fieldname','fe_SSA','data',self.fe_SSA,'format','String'); WriteData(fid,prefix,'object',self,'fieldname','fe_HO' ,'data',self.fe_HO,'format','String'); @@ -182,5 +190,7 @@ writejsdouble(fid,[modelname '.flowequation.isHO'],self.isHO); writejsdouble(fid,[modelname '.flowequation.isFS'],self.isFS); - writejsstring(fid,[modelname '.flowequation.fe_SSA'],self.fe_SSA); + writejsstring(fid,[modelname '.flowequation.isNitscheBC'],self.isNitscheBC); + writejsstring(fid,[modelname '.flowequation.FSNitscheGamma'],self.FSNitscheGamma); + writejsstring(fid,[modelname '.flowequation.fe_SSA'],self.fe_SSA); writejsstring(fid,[modelname '.flowequation.fe_HO'],self.fe_HO); writejsstring(fid,[modelname '.flowequation.fe_FS'],self.fe_FS); Index: /issm/trunk-jpl/src/m/classes/frictionschoof.m =================================================================== --- /issm/trunk-jpl/src/m/classes/frictionschoof.m (revision 25513) +++ /issm/trunk-jpl/src/m/classes/frictionschoof.m (revision 25514) @@ -45,7 +45,7 @@ disp('Schoof sliding law parameters:'); disp(' Schoof''s sliding law reads:'); - disp(' C |u_b|^(m-1) '); + disp(' C^2 |u_b|^(m-1) '); disp(' tau_b = - _____________________________ u_b '); - disp(' (1+(C/(Cmax N))^1/m |u_b| )^m '); + disp(' (1+(C^2/(Cmax N))^1/m |u_b| )^m '); disp(' '); fielddisplay(self,'C','friction coefficient [SI]');