Index: /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp (revision 15435) +++ /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp (revision 15436) @@ -6816,5 +6816,5 @@ /*compute all stiffness matrices for this element*/ Ke1=CreateKMatrixDiagnosticStokesGLSViscous(); - Ke2=CreateKMatrixDiagnosticStokesGLSFriction(); + Ke2=CreateKMatrixDiagnosticStokesFriction(); Ke =new ElementMatrix(Ke1,Ke2); break; @@ -6894,6 +6894,8 @@ ElementMatrix* Penta::CreateKMatrixDiagnosticStokesGLSViscous(void){ + int numdof = NDOF4*NUMVERTICES; + /*Intermediaries */ - int i,approximation; + int i,j,approximation; IssmDouble Jdet,viscosity,stokesreconditioning; IssmDouble xyz_list[NUMVERTICES][3]; @@ -6903,5 +6905,5 @@ IssmDouble D_scalar; IssmDouble D[8][8]={0.0}; - IssmDouble Ke_temp[24][24]={1.0}; //for the six nodes + IssmDouble Ke_temp[24][24]={0.0}; //for the six nodes GaussPenta *gauss=NULL; @@ -6939,4 +6941,6 @@ &B_prime[0][0],8,24,0, &Ke_temp[0][0],1); + + for(i=0;ivalues[i*numdof+j]+=Ke_temp[i][j]; } @@ -6951,76 +6955,4 @@ /*FUNCTION Penta::CreateKMatrixDiagnosticStokesFriction{{{*/ ElementMatrix* Penta::CreateKMatrixDiagnosticStokesFriction(void){ - - /*Constants*/ - const int numdof=NUMVERTICES*NDOF4; - const int numdof2d=NUMVERTICES2D*NDOF4; - - /*Intermediaries */ - int i,j; - int analysis_type,approximation; - IssmDouble alpha2,Jdet2d; - IssmDouble stokesreconditioning,viscosity; - IssmDouble epsilon[6]; /* epsilon=[exx,eyy,ezz,exy,exz,eyz];*/ - IssmDouble xyz_list[NUMVERTICES][3]; - IssmDouble xyz_list_tria[NUMVERTICES2D][3]; - IssmDouble LStokes[2][numdof2d]; - IssmDouble DLStokes[2][2]={0.0}; - IssmDouble Ke_drag_gaussian[numdof2d][numdof2d]; - Friction* friction=NULL; - GaussPenta *gauss=NULL; - - /*If on water or not Stokes, skip stiffness: */ - inputs->GetInputValue(&approximation,ApproximationEnum); - if(IsFloating() || !IsOnBed() || (approximation!=StokesApproximationEnum && approximation!=MacAyealStokesApproximationEnum && approximation!=PattynStokesApproximationEnum)) return NULL; - ElementMatrix* Ke=new ElementMatrix(nodes,NUMVERTICES,this->parameters,StokesApproximationEnum); - - /*Retrieve all inputs and parameters*/ - GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES); - parameters->FindParam(&analysis_type,AnalysisTypeEnum); - parameters->FindParam(&stokesreconditioning,DiagnosticStokesreconditioningEnum); - Input* vx_input=inputs->GetInput(VxEnum); _assert_(vx_input); - Input* vy_input=inputs->GetInput(VyEnum); _assert_(vy_input); - Input* vz_input=inputs->GetInput(VzEnum); _assert_(vz_input); - for(i=0;ibegin();igend();ig++){ - - gauss->GaussPoint(ig); - - GetTriaJacobianDeterminant(&Jdet2d, &xyz_list_tria[0][0],gauss); - GetLStokes(&LStokes[0][0], gauss); - - this->GetStrainRate3d(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input,vz_input); - material->GetViscosity3dStokes(&viscosity,&epsilon[0]); - - friction->GetAlpha2(&alpha2, gauss,VxEnum,VyEnum,VzEnum); - - DLStokes[0][0] = +alpha2*gauss->weight*Jdet2d; //taub_x = -alpha2 vx - DLStokes[1][1] = +alpha2*gauss->weight*Jdet2d; //taub_y = -alpha2 vy - - TripleMultiply( &LStokes[0][0],2,numdof2d,1, - &DLStokes[0][0],2,2,0, - &LStokes[0][0],2,numdof2d,0, - &Ke_drag_gaussian[0][0],0); - - for(i=0;ivalues[i*numdof+j]+=Ke_drag_gaussian[i][j]; - } - - /*DO NOT Transform Coordinate System: this stiffness matrix is already expressed in tangential coordinates*/ - //TransformStiffnessMatrixCoord(Ke,nodes,NUMVERTICES,XYZPEnum); - - /*Clean up and return*/ - delete gauss; - delete friction; - return Ke; -} -/*}}}*/ -/*FUNCTION Penta::CreateKMatrixDiagnosticStokesGLSFriction{{{*/ -ElementMatrix* Penta::CreateKMatrixDiagnosticStokesGLSFriction(void){ /*Constants*/ @@ -7740,8 +7672,27 @@ ElementVector* Penta::CreatePVectorDiagnosticStokes(void){ - /*compute all load vectors for this element*/ - ElementVector* pe1=CreatePVectorDiagnosticStokesViscous(); - ElementVector* pe2=CreatePVectorDiagnosticStokesShelf(); - ElementVector* pe =new ElementVector(pe1,pe2); + int fe_stokes; + ElementVector* pe1; + ElementVector* pe2; + ElementVector* pe; + parameters->FindParam(&fe_stokes,FlowequationFeStokesEnum); + + switch(fe_stokes){ + case 0: + /*compute all stiffness matrices for this element*/ + pe1=CreatePVectorDiagnosticStokesViscous(); + pe2=CreatePVectorDiagnosticStokesShelf(); + pe =new ElementVector(pe1,pe2); + break; + case 1: + /*compute all stiffness matrices for this element*/ + pe1=CreatePVectorDiagnosticStokesGLSViscous(); + pe2=CreatePVectorDiagnosticStokesShelf(); + pe =new ElementVector(pe1,pe2); + break; + default: + _error_("Finite element" << fe_stokes << " not supported yet"); + } + /*clean-up and return*/ @@ -7830,6 +7781,123 @@ } + if(IsOnBed() || IsOnSurface()){ + IssmDouble xyz_list_tria[NUMVERTICES2D][3]; + IssmDouble pi=3.141592653589793; + IssmDouble x,y,z; + IssmDouble basis[6]; //for the six nodes of the penta + GaussPenta *gauss=NULL; + Input* x_input=inputs->GetInput(MeshXEnum); _assert_(x_input); + Input* y_input=inputs->GetInput(MeshYEnum); _assert_(y_input); + Input* z_input=inputs->GetInput(MeshZEnum); _assert_(z_input); + + + if(IsOnBed()){ + for(i=0;ibegin();igend();ig++){ + + gauss->GaussPoint(ig); + + GetTriaJacobianDeterminant(&Jdet, &xyz_list_tria[0][0], gauss); + GetNodalFunctionsP1(basis, gauss); + x_input->GetInputValue(&x, gauss); + y_input->GetInputValue(&y, gauss); + z_input->GetInputValue(&z, gauss); + + + forcex=-((cos(2*pi*z)-1)*2*pi*sin(2*pi*y)*cos(2*pi*x)-2*(cos(2*pi*x)-1)*2*pi*sin(2*pi*y)*cos(2*pi*z))*1; + forcey=-((cos(2*pi*z)-1)*2*pi*sin(2*pi*x)*cos(2*pi*y)+(cos(2*pi*y)-1)*2*pi*sin(2*pi*x)*cos(2*pi*z))*1; + forcez=-(-2*pi*2*sin(2*pi*x)*sin(2*pi*y)*sin(2*pi*z)+sin(2*pi*x)*sin(2*pi*y)*sin(2*pi*z)); + + if(IsOnBed()){ + for(i=0;i<3;i++){ + Pe_gaussian[i*NDOF4+0]+=forcex*Jdet*gauss->weight*basis[i]*-1.; + Pe_gaussian[i*NDOF4+1]+=forcey*Jdet*gauss->weight*basis[i]*-1.; + Pe_gaussian[i*NDOF4+2]+=forcez*Jdet*gauss->weight*basis[i]*-1.; + } + } + if(IsOnSurface()){ + for(i=3;i<6;i++){ + Pe_gaussian[i*NDOF4+0]+=forcex*Jdet*gauss->weight*basis[i]*1.; + Pe_gaussian[i*NDOF4+1]+=forcey*Jdet*gauss->weight*basis[i]*1.; + Pe_gaussian[i*NDOF4+2]+=forcez*Jdet*gauss->weight*basis[i]*1.; + } + } + } + } + /*Condensation*/ ReduceVectorStokes(pe->values, &Ke_temp[0][0], &Pe_gaussian[0]); + + /*Transform coordinate system*/ + TransformLoadVectorCoord(pe,nodes,NUMVERTICES,XYZPEnum); + + /*Clean up and return*/ + delete gauss; + return pe; +} +/*}}}*/ +/*FUNCTION Penta::CreatePVectorDiagnosticStokesGLSViscous {{{*/ +ElementVector* Penta::CreatePVectorDiagnosticStokesGLSViscous(void){ + + /*Constants*/ + const int numdof=NDOF4*NUMVERTICES; + + /*Intermediaries*/ + int i,j; + int approximation; + IssmDouble Jdet,gravity,rho_ice; + IssmDouble forcex,forcey,forcez; + IssmDouble xyz_list[NUMVERTICES][3]; + IssmDouble epsilon[6]; /* epsilon=[exx,eyy,ezz,exy,exz,eyz];*/ + IssmDouble l1l6[6]; //for the six nodes and the bubble + IssmDouble Pe_gaussian[numdof]={0.0}; //for the six nodes and the bubble + GaussPenta *gauss=NULL; + + /*Initialize Element vector and return if necessary*/ + inputs->GetInputValue(&approximation,ApproximationEnum); + if(approximation!=StokesApproximationEnum && approximation!=MacAyealStokesApproximationEnum && approximation!=PattynStokesApproximationEnum) return NULL; + ElementVector* pe=new ElementVector(nodes,NUMVERTICES,this->parameters,StokesApproximationEnum); + + /*Retrieve all inputs and parameters*/ + rho_ice=matpar->GetRhoIce(); + gravity=matpar->GetG(); + GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES); + Input* vx_input=inputs->GetInput(VxEnum); _assert_(vx_input); + Input* vy_input=inputs->GetInput(VyEnum); _assert_(vy_input); + Input* vz_input=inputs->GetInput(VzEnum); _assert_(vz_input); + Input* loadingforcex_input=inputs->GetInput(LoadingforceXEnum); _assert_(loadingforcex_input); + Input* loadingforcey_input=inputs->GetInput(LoadingforceYEnum); _assert_(loadingforcey_input); + Input* loadingforcez_input=inputs->GetInput(LoadingforceZEnum); _assert_(loadingforcez_input); + + /* Start looping on the number of gaussian points: */ + gauss=new GaussPenta(5,5); + for(int ig=gauss->begin();igend();ig++){ + + gauss->GaussPoint(ig); + + GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss); + GetNodalFunctionsP1(&l1l6[0], gauss); + + loadingforcex_input->GetInputValue(&forcex, gauss); + loadingforcey_input->GetInputValue(&forcey, gauss); + loadingforcez_input->GetInputValue(&forcez, gauss); + + for(i=0;iweight*l1l6[i]; + Pe_gaussian[i*NDOF4+0]+=forcex*Jdet*gauss->weight*l1l6[i]; + Pe_gaussian[i*NDOF4+1]+=forcey*Jdet*gauss->weight*l1l6[i]; + Pe_gaussian[i*NDOF4+2]+=forcez*Jdet*gauss->weight*l1l6[i]; + } + } + + for(i=0;ivalues[i*NDOF4+j]+=Pe_gaussian[i*NDOF4+j]; /*Transform coordinate system*/ Index: /issm/trunk-jpl/src/c/classes/Elements/Penta.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Penta.h (revision 15435) +++ /issm/trunk-jpl/src/c/classes/Elements/Penta.h (revision 15436) @@ -287,4 +287,5 @@ ElementVector* CreatePVectorDiagnosticStokes(void); ElementVector* CreatePVectorDiagnosticStokesViscous(void); + ElementVector* CreatePVectorDiagnosticStokesGLSViscous(void); ElementVector* CreatePVectorDiagnosticStokesShelf(void); ElementVector* CreatePVectorDiagnosticVert(void); Index: /issm/trunk-jpl/src/c/classes/Elements/PentaRef.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/PentaRef.cpp (revision 15435) +++ /issm/trunk-jpl/src/c/classes/Elements/PentaRef.cpp (revision 15436) @@ -275,28 +275,28 @@ /*Build B: */ for (i=0;i