Index: /issm/trunk/src/c/objects/Elements/Penta.cpp =================================================================== --- /issm/trunk/src/c/objects/Elements/Penta.cpp (revision 5852) +++ /issm/trunk/src/c/objects/Elements/Penta.cpp (revision 5853) @@ -2214,5 +2214,7 @@ break; case StokesApproximationEnum: - CreateKMatrixDiagnosticStokes( Kgg); + Ke=CreateKMatrixDiagnosticStokes(); + if(Ke) Ke->AddToGlobal(Kgg,NULL,NULL); + delete Ke; break; case HutterApproximationEnum: @@ -2233,5 +2235,7 @@ if(Ke) Ke->AddToGlobal(Kgg,NULL,NULL); delete Ke; - CreateKMatrixDiagnosticStokes( Kgg); + Ke=CreateKMatrixDiagnosticStokes(); + if(Ke) Ke->AddToGlobal(Kgg,NULL,NULL); + delete Ke; CreateKMatrixCouplingPattynStokes( Kgg); break; @@ -2556,19 +2560,104 @@ } /*}}}*/ -/*FUNCTION Penta::CreateKMatrixDiagnosticStokes {{{1*/ -void Penta::CreateKMatrixDiagnosticStokes( Mat Kgg){ - +/*FUNCTION Penta::CreateKMatrixDiagnosticStokes{{{1*/ +ElementMatrix* Penta::CreateKMatrixDiagnosticStokes(void){ + + /*compute all stiffness matrices for this element*/ + ElementMatrix* Ke1=CreateKMatrixDiagnosticStokesViscous(); + ElementMatrix* Ke2=CreateKMatrixDiagnosticStokesFriction(); + ElementMatrix* Ke =new ElementMatrix(Ke1,Ke2); + + /*clean-up and return*/ + delete Ke1; + delete Ke2; + return Ke; +} +/*}}}*/ +/*FUNCTION Penta::CreateKMatrixDiagnosticStokesViscous {{{1*/ +ElementMatrix* Penta::CreateKMatrixDiagnosticStokesViscous(void){ + + const int numdof=NUMVERTICES*NDOF4; int i,j; - const int numdof=NUMVERTICES*NDOF4; - int* doflist=NULL; - - const int numdof2d=NUMVERTICES2D*NDOF4; - - /*Grid data: */ + int ig; double xyz_list[NUMVERTICES][3]; double xyz_list_tria[NUMVERTICES2D][3]; double bed_normal[3]; - - /*matrices: */ + double Ke_temp[27][27]={0.0}; //for the six nodes and the bubble + double Ke_reduced[numdof][numdof]; //for the six nodes only + double Ke_gaussian[27][27]; + double B[8][27]; + double B_prime[8][27]; + double Jdet; + double D[8][8]={0.0}; + double D_scalar; + double DLStokes[14][14]={0.0}; + GaussPenta *gauss=NULL; + double epsilon[6]; /* epsilon=[exx,eyy,ezz,exy,exz,eyz];*/ + double viscosity; + double alpha2_gauss; + Friction* friction=NULL; + double stokesreconditioning; + int analysis_type; + int approximation; + + /*If on water or not Stokes, skip stiffness: */ + inputs->GetParameterValue(&approximation,ApproximationEnum); + if(IsOnWater() || approximation!=StokesApproximationEnum) return NULL; + ElementMatrix* Ke=this->NewElementMatrix(StokesApproximationEnum); + + /*Retrieve all inputs and parameters*/ + GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES); + parameters->FindParam(&analysis_type,AnalysisTypeEnum); + parameters->FindParam(&stokesreconditioning,StokesReconditioningEnum); + Input* vx_input=inputs->GetInput(VxEnum); ISSMASSERT(vx_input); + Input* vy_input=inputs->GetInput(VyEnum); ISSMASSERT(vy_input); + Input* vz_input=inputs->GetInput(VzEnum); ISSMASSERT(vz_input); + + /* Start looping on the number of gaussian points: */ + gauss=new GaussPenta(5,5); + for (ig=gauss->begin();igend();ig++){ + + gauss->GaussPoint(ig); + + this->GetStrainRate3d(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input,vz_input); + matice->GetViscosity3dStokes(&viscosity,&epsilon[0]); + + GetBStokes(&B[0][0],&xyz_list[0][0],gauss); + GetBprimeStokes(&B_prime[0][0],&xyz_list[0][0],gauss); + + GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss); + + /* Build the D matrix: we plug the gaussian weight, the thickness, the viscosity, and the jacobian determinant + * onto this scalar matrix, so that we win some computational time: */ + D_scalar=gauss->weight*Jdet; + for (i=0;i<6;i++) D[i][i]=D_scalar*2*viscosity; + for (i=6;i<8;i++) D[i][i]=-D_scalar*stokesreconditioning; + + TripleMultiply( &B[0][0],8,27,1, + &D[0][0],8,8,0, + &B_prime[0][0],8,27,0, + &Ke_gaussian[0][0],0); + + for(i=0;i<27;i++) for(j=0;j<27;j++) Ke_temp[i][j]+=Ke_gaussian[i][j]; + } + + /*Condensation*/ + ReduceMatrixStokes(Ke->values, &Ke_temp[0][0]); + + /*Clean up and return*/ + delete gauss; + return Ke; +} +/*}}}*/ +/*FUNCTION Penta::CreateKMatrixDiagnosticStokesFriction {{{1*/ +ElementMatrix* Penta::CreateKMatrixDiagnosticStokesFriction(void){ + + const int numdof=NUMVERTICES*NDOF4; + const int numdof2d=NUMVERTICES2D*NDOF4; + int i,j; + int ig; + double xyz_list[NUMVERTICES][3]; + double xyz_list_tria[NUMVERTICES2D][3]; + double bed_normal[3]; double Ke_temp[27][27]={0.0}; //for the six nodes and the bubble double Ke_reduced[numdof][numdof]; //for the six nodes only @@ -2584,7 +2673,4 @@ double D_scalar; double DLStokes[14][14]={0.0}; - - /* gaussian points: */ - int ig; GaussPenta *gauss=NULL; double epsilon[6]; /* epsilon=[exx,eyy,ezz,exy,exz,eyz];*/ @@ -2592,139 +2678,70 @@ double alpha2_gauss; Friction* friction=NULL; - - /*parameters: */ double stokesreconditioning; int analysis_type; int approximation; - /*retrive parameters: */ + /*If on water or not Stokes, skip stiffness: */ + inputs->GetParameterValue(&approximation,ApproximationEnum); + if(IsOnWater() || IsOnShelf() || !IsOnBed() || approximation!=StokesApproximationEnum) return NULL; + ElementMatrix* Ke=this->NewElementMatrix(StokesApproximationEnum); + + /*Retrieve all inputs and parameters*/ + GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES); parameters->FindParam(&analysis_type,AnalysisTypeEnum); - this->parameters->FindParam(&stokesreconditioning,StokesReconditioningEnum); - - /*retrieve inputs :*/ - inputs->GetParameterValue(&approximation,ApproximationEnum); - - /*If on water or not Stokes, skip stiffness: */ - if(IsOnWater() || approximation!=StokesApproximationEnum) return; - - /* Get node coordinates and dof list: */ - GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES); - GetDofList(&doflist,StokesApproximationEnum,GsetEnum); - - /*Retrieve all inputs we will be needing: */ + parameters->FindParam(&stokesreconditioning,StokesReconditioningEnum); Input* vx_input=inputs->GetInput(VxEnum); ISSMASSERT(vx_input); Input* vy_input=inputs->GetInput(VyEnum); ISSMASSERT(vy_input); Input* vz_input=inputs->GetInput(VzEnum); ISSMASSERT(vz_input); + for(i=0;ibegin();igend();ig++){ gauss->GaussPoint(ig); - /*Compute strain rate: */ + GetLStokes(&LStokes[0][0], gauss); + GetLprimeStokes(&LprimeStokes[0][0], &xyz_list[0][0], gauss); + this->GetStrainRate3d(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input,vz_input); - - /*Get viscosity: */ matice->GetViscosity3dStokes(&viscosity,&epsilon[0]); - /*Get B and Bprime matrices: */ - GetBStokes(&B[0][0],&xyz_list[0][0],gauss); - GetBprimeStokes(&B_prime[0][0],&xyz_list[0][0],gauss); - - /* Get Jacobian determinant: */ - GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss); - - /* Build the D matrix: we plug the gaussian weight, the thickness, the viscosity, and the jacobian determinant - * onto this scalar matrix, so that we win some computational time: */ - D_scalar=gauss->weight*Jdet; - for (i=0;i<6;i++) D[i][i]=D_scalar*2*viscosity; - for (i=6;i<8;i++) D[i][i]=-D_scalar*stokesreconditioning; - - /* Do the triple product tB*D*Bprime: */ - TripleMultiply( &B[0][0],8,27,1, - &D[0][0],8,8,0, - &B_prime[0][0],8,27,0, - &Ke_gaussian[0][0],0); - - /*Add Ke_gaussian and Ke_gaussian to terms in pKe. Watch out for column orientation from matlab: */ - for(i=0;i<27;i++) for(j=0;j<27;j++) Ke_temp[i][j]+=Ke_gaussian[i][j]; - } - delete gauss; //gauss of previous loop - - if(IsOnBed() && !IsOnShelf()){ - - /*build friction object, used later on: */ - friction=new Friction("3d",inputs,matpar,analysis_type); - - for(i=0;ibegin();igend();ig++){ - - gauss->GaussPoint(ig); - - /*Get the Jacobian determinant */ - GetTriaJacobianDeterminant(&Jdet2d, &xyz_list_tria[0][0],gauss); - - /*Get L matrix if viscous basal drag present: */ - GetLStokes(&LStokes[0][0], gauss); - GetLprimeStokes(&LprimeStokes[0][0], &xyz_list[0][0], gauss); - - /*Compute strain rate: */ - this->GetStrainRate3d(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input,vz_input); - - /*Get viscosity at last iteration: */ - matice->GetViscosity3dStokes(&viscosity,&epsilon[0]); - - /*Get normal vecyor to the bed */ - BedNormal(&bed_normal[0],xyz_list_tria); - - /*Calculate DL on gauss point */ - friction->GetAlpha2(&alpha2_gauss, gauss,VxEnum,VyEnum,VzEnum); - - DLStokes[0][0]=alpha2_gauss*gauss->weight*Jdet2d; - DLStokes[1][1]=alpha2_gauss*gauss->weight*Jdet2d; - DLStokes[2][2]=-alpha2_gauss*gauss->weight*Jdet2d*bed_normal[0]*bed_normal[2]; - DLStokes[3][3]=-alpha2_gauss*gauss->weight*Jdet2d*bed_normal[1]*bed_normal[2]; - DLStokes[4][4]=-alpha2_gauss*gauss->weight*Jdet2d*bed_normal[0]*bed_normal[2]; - DLStokes[5][5]=-alpha2_gauss*gauss->weight*Jdet2d*bed_normal[1]*bed_normal[2]; - DLStokes[6][6]=-2*viscosity*gauss->weight*Jdet2d*bed_normal[0]; - DLStokes[7][7]=-2*viscosity*gauss->weight*Jdet2d*bed_normal[1]; - DLStokes[8][8]=-2*viscosity*gauss->weight*Jdet2d*bed_normal[2]; - DLStokes[9][8]=-2*viscosity*gauss->weight*Jdet2d*bed_normal[0]/2.0; - DLStokes[10][10]=-2*viscosity*gauss->weight*Jdet2d*bed_normal[1]/2.0; - DLStokes[11][11]=stokesreconditioning*gauss->weight*Jdet2d*bed_normal[0]; - DLStokes[12][12]=stokesreconditioning*gauss->weight*Jdet2d*bed_normal[1]; - DLStokes[13][13]=stokesreconditioning*gauss->weight*Jdet2d*bed_normal[2]; - - /* Do the triple product tL*D*L: */ - TripleMultiply( &LStokes[0][0],14,numdof2d,1, - &DLStokes[0][0],14,14,0, - &LprimeStokes[0][0],14,numdof2d,0, - &Ke_drag_gaussian[0][0],0); - - for(i=0;iGetAlpha2(&alpha2_gauss, gauss,VxEnum,VyEnum,VzEnum); + + DLStokes[0][0]=alpha2_gauss*gauss->weight*Jdet2d; + DLStokes[1][1]=alpha2_gauss*gauss->weight*Jdet2d; + DLStokes[2][2]=-alpha2_gauss*gauss->weight*Jdet2d*bed_normal[0]*bed_normal[2]; + DLStokes[3][3]=-alpha2_gauss*gauss->weight*Jdet2d*bed_normal[1]*bed_normal[2]; + DLStokes[4][4]=-alpha2_gauss*gauss->weight*Jdet2d*bed_normal[0]*bed_normal[2]; + DLStokes[5][5]=-alpha2_gauss*gauss->weight*Jdet2d*bed_normal[1]*bed_normal[2]; + DLStokes[6][6]=-2*viscosity*gauss->weight*Jdet2d*bed_normal[0]; + DLStokes[7][7]=-2*viscosity*gauss->weight*Jdet2d*bed_normal[1]; + DLStokes[8][8]=-2*viscosity*gauss->weight*Jdet2d*bed_normal[2]; + DLStokes[9][8]=-2*viscosity*gauss->weight*Jdet2d*bed_normal[0]/2.0; + DLStokes[10][10]=-2*viscosity*gauss->weight*Jdet2d*bed_normal[1]/2.0; + DLStokes[11][11]=stokesreconditioning*gauss->weight*Jdet2d*bed_normal[0]; + DLStokes[12][12]=stokesreconditioning*gauss->weight*Jdet2d*bed_normal[1]; + DLStokes[13][13]=stokesreconditioning*gauss->weight*Jdet2d*bed_normal[2]; + + TripleMultiply( &LStokes[0][0],14,numdof2d,1, + &DLStokes[0][0],14,14,0, + &LprimeStokes[0][0],14,numdof2d,0, + &Ke_drag_gaussian[0][0],0); + + for(i=0;ivalues[i*numdof+j]+=Ke_temp[i][j]; + + /*Clean up and return*/ + delete gauss; + delete friction; + return Ke; } /*}}}*/ Index: /issm/trunk/src/c/objects/Elements/Penta.h =================================================================== --- /issm/trunk/src/c/objects/Elements/Penta.h (revision 5852) +++ /issm/trunk/src/c/objects/Elements/Penta.h (revision 5853) @@ -136,5 +136,7 @@ ElementMatrix* CreateKMatrixDiagnosticPattynViscous(void); ElementMatrix* CreateKMatrixDiagnosticPattynFriction(void); - void CreateKMatrixDiagnosticStokes( Mat Kgg); + ElementMatrix* CreateKMatrixDiagnosticStokes(void); + ElementMatrix* CreateKMatrixDiagnosticStokesViscous(void); + ElementMatrix* CreateKMatrixDiagnosticStokesFriction(void); void CreateKMatrixDiagnosticVert( Mat Kgg); void CreateKMatrixMelting(Mat Kggg);