Index: /issm/trunk/src/c/objects/Elements/Penta.cpp =================================================================== --- /issm/trunk/src/c/objects/Elements/Penta.cpp (revision 7069) +++ /issm/trunk/src/c/objects/Elements/Penta.cpp (revision 7070) @@ -704,5 +704,5 @@ /*compute all stiffness matrices for this element*/ ElementMatrix* Ke1=CreateKMatrixCouplingMacAyealStokesViscous(); - ElementMatrix* Ke2=CreateKMatrixCouplingMacAyealPattynFriction(); + ElementMatrix* Ke2=CreateKMatrixCouplingMacAyealStokesFriction(); ElementMatrix* Ke =new ElementMatrix(Ke1,Ke2); @@ -804,14 +804,103 @@ } /*}}}*/ -/*FUNCTION Penta::CreateKMatrixCouplingMacAyealStokesFriction{{{1*/ +/*FUNCTION Penta::CreateKMatrixCouplingMacAyealStokesFriction {{{1*/ ElementMatrix* Penta::CreateKMatrixCouplingMacAyealStokesFriction(void){ - /*Initialize Element matrix and return if necessary*/ + /*Constants*/ + const int numdof=NUMVERTICES*NDOF4; + const int numdofm=NUMVERTICES*NDOF2; + const int numdof2d=NUMVERTICES2D*NDOF4; + const int numdof2dm=NUMVERTICES2D*NDOF2; + const int numdoftot=numdof+numdofm; + + /*Intermediaries */ + int i,j,ig; + int analysis_type,approximation; + double stokesreconditioning; + double viscosity,alpha2_gauss,Jdet2d; + double bed_normal[3]; + double epsilon[6]; /* epsilon=[exx,eyy,ezz,exy,exz,eyz];*/ + double xyz_list[NUMVERTICES][3]; + double xyz_list_tria[NUMVERTICES2D][3]; + double LMacAyealStokes[8][numdof2dm]; + double LprimeMacAyealStokes[8][numdof2d]; + double DLMacAyealStokes[8][8]={0.0}; + double LStokesMacAyeal[4][numdof2d]; + double LprimeStokesMacAyeal[4][numdof2dm]; + double DLStokesMacAyeal[4][4]={0.0}; + double Ke_drag_gaussian[numdof2dm][numdof2d]; + double Ke_drag_gaussian2[numdof2d][numdof2dm]; + Friction* friction=NULL; + GaussPenta *gauss=NULL; + + /*If on water or not Stokes, skip stiffness: */ + inputs->GetParameterValue(&approximation,ApproximationEnum); if(IsOnWater() || IsOnShelf() || !IsOnBed()) return NULL; - - Tria* tria=(Tria*)SpawnTria(0,1,2); //grids 0, 1 and 2 make the new tria. - ElementMatrix* Ke=tria->CreateKMatrixCouplingMacAyealPattynFriction(); - delete tria->matice; delete tria; - + ElementMatrix* Ke1=new ElementMatrix(this->nodes,NUMVERTICES,this->parameters,MacAyealApproximationEnum); + ElementMatrix* Ke2=new ElementMatrix(this->nodes,NUMVERTICES,this->parameters,StokesApproximationEnum); + ElementMatrix* Ke=new ElementMatrix(Ke1,Ke2); + delete Ke1; delete Ke2; + + /*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); _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); + GetLMacAyealStokes(&LMacAyealStokes[0][0], gauss); + GetLprimeMacAyealStokes(&LprimeMacAyealStokes[0][0], &xyz_list[0][0], gauss); + GetLStokesMacAyeal(&LStokesMacAyeal[0][0], gauss); + GetLprimeStokesMacAyeal(&LprimeStokesMacAyeal[0][0], &xyz_list[0][0], gauss); + + this->GetStrainRate3d(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input,vz_input); + matice->GetViscosity3dStokes(&viscosity,&epsilon[0]); + + BedNormal(&bed_normal[0],xyz_list_tria); + friction->GetAlpha2(&alpha2_gauss, gauss,VxEnum,VyEnum,VzEnum); + + DLMacAyealStokes[0][0]=alpha2_gauss*gauss->weight*Jdet2d; + DLMacAyealStokes[1][1]=alpha2_gauss*gauss->weight*Jdet2d; + DLMacAyealStokes[2][2]=-alpha2_gauss*gauss->weight*Jdet2d*bed_normal[0]*bed_normal[2]; + DLMacAyealStokes[3][3]=-alpha2_gauss*gauss->weight*Jdet2d*bed_normal[1]*bed_normal[2]; + DLMacAyealStokes[4][4]=-2*viscosity*gauss->weight*Jdet2d*bed_normal[0]; + DLMacAyealStokes[5][5]=-2*viscosity*gauss->weight*Jdet2d*bed_normal[1]; + DLMacAyealStokes[6][6]=stokesreconditioning*gauss->weight*Jdet2d*bed_normal[0]; + DLMacAyealStokes[7][7]=stokesreconditioning*gauss->weight*Jdet2d*bed_normal[1]; + + DLStokesMacAyeal[0][0]=alpha2_gauss*gauss->weight*Jdet2d; + DLStokesMacAyeal[1][1]=alpha2_gauss*gauss->weight*Jdet2d; + DLStokesMacAyeal[2][2]=-alpha2_gauss*gauss->weight*Jdet2d*bed_normal[0]*bed_normal[2]; + DLStokesMacAyeal[3][3]=-alpha2_gauss*gauss->weight*Jdet2d*bed_normal[1]*bed_normal[2]; + + TripleMultiply( &LMacAyealStokes[0][0],8,numdof2dm,1, + &DLMacAyealStokes[0][0],8,8,0, + &LprimeMacAyealStokes[0][0],8,numdof2d,0, + &Ke_drag_gaussian[0][0],0); + + TripleMultiply( &LStokesMacAyeal[0][0],4,numdof2d,1, + &DLStokesMacAyeal[0][0],4,4,0, + &LprimeStokesMacAyeal[0][0],4,numdof2dm,0, + &Ke_drag_gaussian2[0][0],0); + + for(i=0;ivalues[i*numdoftot+j+numdofm]+=Ke_drag_gaussian[i][j]; + for(i=0;ivalues[(i+numdofm)*numdoftot+j]+=Ke_drag_gaussian2[i][j]; + } + + /*Clean up and return*/ + delete gauss; + delete friction; return Ke; } Index: /issm/trunk/src/c/objects/Elements/PentaRef.cpp =================================================================== --- /issm/trunk/src/c/objects/Elements/PentaRef.cpp (revision 7069) +++ /issm/trunk/src/c/objects/Elements/PentaRef.cpp (revision 7070) @@ -713,4 +713,206 @@ } /*}}}*/ +/*FUNCTION PentaRef::GetLMacAyealStokes {{{1*/ +void PentaRef::GetLMacAyealStokes(double* LStokes, GaussPenta* gauss){ + /* + * Compute L matrix. L=[L1 L2 L3] where Li is square and of size numdof. + * For grid i, Li can be expressed in the actual coordinate system + * by: + * Li=[ h 0 ] + * [ 0 h ] + * [ h 0 ] + * [ 0 h ] + * [ h 0 ] + * [ 0 h ] + * [ h 0 ] + * [ 0 h ] + * where h is the interpolation function for grid i. + */ + + int i; + int num_dof=2; + + double l1l2l3[NUMNODESP1_2d]; + + + /*Get l1l2l3 in actual coordinate system: */ + l1l2l3[0]=gauss->coord1*(1-gauss->coord4)/2.0; + l1l2l3[1]=gauss->coord2*(1-gauss->coord4)/2.0; + l1l2l3[2]=gauss->coord3*(1-gauss->coord4)/2.0; + + /*Build LStokes: */ + for (i=0;i<3;i++){ + *(LStokes+num_dof*NUMNODESP1_2d*0+num_dof*i)=l1l2l3[i]; //LStokes[0][NDOF2*i]=dh1dh3[0][i]; + *(LStokes+num_dof*NUMNODESP1_2d*0+num_dof*i+1)=0; + *(LStokes+num_dof*NUMNODESP1_2d*1+num_dof*i)=0; + *(LStokes+num_dof*NUMNODESP1_2d*1+num_dof*i+1)=l1l2l3[i]; + *(LStokes+num_dof*NUMNODESP1_2d*2+num_dof*i)=l1l2l3[i]; + *(LStokes+num_dof*NUMNODESP1_2d*2+num_dof*i+1)=0; + *(LStokes+num_dof*NUMNODESP1_2d*3+num_dof*i)=0; + *(LStokes+num_dof*NUMNODESP1_2d*3+num_dof*i+1)=l1l2l3[i]; + *(LStokes+num_dof*NUMNODESP1_2d*4+num_dof*i)=l1l2l3[i]; + *(LStokes+num_dof*NUMNODESP1_2d*4+num_dof*i+1)=0; + *(LStokes+num_dof*NUMNODESP1_2d*5+num_dof*i)=0; + *(LStokes+num_dof*NUMNODESP1_2d*5+num_dof*i+1)=l1l2l3[i]; + *(LStokes+num_dof*NUMNODESP1_2d*6+num_dof*i)=l1l2l3[i]; + *(LStokes+num_dof*NUMNODESP1_2d*6+num_dof*i+1)=0; + *(LStokes+num_dof*NUMNODESP1_2d*7+num_dof*i)=0; + *(LStokes+num_dof*NUMNODESP1_2d*7+num_dof*i+1)=l1l2l3[i]; + + } +} +/*}}}*/ +/*FUNCTION PentaRef::GetLprimeMacAyealStokes {{{1*/ +void PentaRef::GetLprimeMacAyealStokes(double* LprimeStokes, double* xyz_list, GaussPenta* gauss){ + + /* + * Compute Lprime matrix. Lprime=[Lp1 Lp2 Lp3] where Lpi is square and of size numdof. + * For grid i, Lpi can be expressed in the actual coordinate system + * by: + * Lpi=[ h 0 0 0] + * [ 0 h 0 0] + * [ 0 0 h 0] + * [ 0 0 h 0] + * [ 0 0 dh/dz 0] + * [ 0 0 dh/dz 0] + * [ 0 0 0 h] + * [ 0 0 0 h] + * where h is the interpolation function for grid i. + */ + int i; + int num_dof=4; + + double l1l2l3[NUMNODESP1_2d]; + double dh1dh6[3][NUMNODESP1]; + + /*Get l1l2l3 in actual coordinate system: */ + l1l2l3[0]=gauss->coord1*(1-gauss->coord4)/2.0; + l1l2l3[1]=gauss->coord2*(1-gauss->coord4)/2.0; + l1l2l3[2]=gauss->coord3*(1-gauss->coord4)/2.0; + + GetNodalFunctionsP1Derivatives(&dh1dh6[0][0],xyz_list,gauss); + + /*Build LprimeStokes: */ + for (i=0;i<3;i++){ + *(LprimeStokes+num_dof*NUMNODESP1_2d*0+num_dof*i)=l1l2l3[i]; //LprimeStokes[0][NDOF2*i]=dh1dh3[0][i]; + *(LprimeStokes+num_dof*NUMNODESP1_2d*0+num_dof*i+1)=0; + *(LprimeStokes+num_dof*NUMNODESP1_2d*0+num_dof*i+2)=0; + *(LprimeStokes+num_dof*NUMNODESP1_2d*0+num_dof*i+3)=0; + *(LprimeStokes+num_dof*NUMNODESP1_2d*1+num_dof*i)=0; + *(LprimeStokes+num_dof*NUMNODESP1_2d*1+num_dof*i+1)=l1l2l3[i]; + *(LprimeStokes+num_dof*NUMNODESP1_2d*1+num_dof*i+2)=0; + *(LprimeStokes+num_dof*NUMNODESP1_2d*1+num_dof*i+3)=0; + *(LprimeStokes+num_dof*NUMNODESP1_2d*2+num_dof*i)=0; + *(LprimeStokes+num_dof*NUMNODESP1_2d*2+num_dof*i+1)=0; + *(LprimeStokes+num_dof*NUMNODESP1_2d*2+num_dof*i+2)=l1l2l3[i]; + *(LprimeStokes+num_dof*NUMNODESP1_2d*2+num_dof*i+3)=0; + *(LprimeStokes+num_dof*NUMNODESP1_2d*3+num_dof*i)=0; + *(LprimeStokes+num_dof*NUMNODESP1_2d*3+num_dof*i+1)=0; + *(LprimeStokes+num_dof*NUMNODESP1_2d*3+num_dof*i+2)=l1l2l3[i]; + *(LprimeStokes+num_dof*NUMNODESP1_2d*3+num_dof*i+3)=0; + *(LprimeStokes+num_dof*NUMNODESP1_2d*4+num_dof*i)=0; + *(LprimeStokes+num_dof*NUMNODESP1_2d*4+num_dof*i+1)=0; + *(LprimeStokes+num_dof*NUMNODESP1_2d*4+num_dof*i+2)=dh1dh6[2][i]; + *(LprimeStokes+num_dof*NUMNODESP1_2d*4+num_dof*i+3)=0; + *(LprimeStokes+num_dof*NUMNODESP1_2d*5+num_dof*i)=0; + *(LprimeStokes+num_dof*NUMNODESP1_2d*5+num_dof*i+1)=0; + *(LprimeStokes+num_dof*NUMNODESP1_2d*5+num_dof*i+2)=dh1dh6[2][i]; + *(LprimeStokes+num_dof*NUMNODESP1_2d*5+num_dof*i+3)=0; + *(LprimeStokes+num_dof*NUMNODESP1_2d*6+num_dof*i)=0; + *(LprimeStokes+num_dof*NUMNODESP1_2d*6+num_dof*i+1)=0; + *(LprimeStokes+num_dof*NUMNODESP1_2d*6+num_dof*i+2)=0; + *(LprimeStokes+num_dof*NUMNODESP1_2d*6+num_dof*i+3)=l1l2l3[i]; + *(LprimeStokes+num_dof*NUMNODESP1_2d*7+num_dof*i)=0; + *(LprimeStokes+num_dof*NUMNODESP1_2d*7+num_dof*i+1)=0; + *(LprimeStokes+num_dof*NUMNODESP1_2d*7+num_dof*i+2)=0; + *(LprimeStokes+num_dof*NUMNODESP1_2d*7+num_dof*i+3)=l1l2l3[i]; + } +} +/*}}}*/ +/*FUNCTION PentaRef::GetLStokesMacAyeal {{{1*/ +void PentaRef::GetLStokesMacAyeal(double* LStokes, GaussPenta* gauss){ + /* + * Compute L matrix. L=[L1 L2 L3] where Li is square and of size numdof. + * For grid i, Li can be expressed in the actual coordinate system + * by: + * Li=[ h 0 0 0] + * [ 0 h 0 0] + * [ 0 0 h 0] + * [ 0 0 h 0] + * where h is the interpolation function for grid i. + */ + + int i; + int num_dof=4; + + double l1l2l3[NUMNODESP1_2d]; + + + /*Get l1l2l3 in actual coordinate system: */ + l1l2l3[0]=gauss->coord1*(1-gauss->coord4)/2.0; + l1l2l3[1]=gauss->coord2*(1-gauss->coord4)/2.0; + l1l2l3[2]=gauss->coord3*(1-gauss->coord4)/2.0; + + /*Build LStokes: */ + for (i=0;i<3;i++){ + *(LStokes+num_dof*NUMNODESP1_2d*0+num_dof*i)=l1l2l3[i]; //LStokes[0][NDOF2*i]=dh1dh3[0][i]; + *(LStokes+num_dof*NUMNODESP1_2d*0+num_dof*i+1)=0; + *(LStokes+num_dof*NUMNODESP1_2d*0+num_dof*i+2)=0; + *(LStokes+num_dof*NUMNODESP1_2d*0+num_dof*i+3)=0; + *(LStokes+num_dof*NUMNODESP1_2d*1+num_dof*i)=0; + *(LStokes+num_dof*NUMNODESP1_2d*1+num_dof*i+1)=l1l2l3[i]; + *(LStokes+num_dof*NUMNODESP1_2d*1+num_dof*i+2)=0; + *(LStokes+num_dof*NUMNODESP1_2d*1+num_dof*i+3)=0; + *(LStokes+num_dof*NUMNODESP1_2d*2+num_dof*i)=0; + *(LStokes+num_dof*NUMNODESP1_2d*2+num_dof*i+1)=0; + *(LStokes+num_dof*NUMNODESP1_2d*2+num_dof*i+2)=l1l2l3[i]; + *(LStokes+num_dof*NUMNODESP1_2d*2+num_dof*i+3)=0; + *(LStokes+num_dof*NUMNODESP1_2d*3+num_dof*i)=0; + *(LStokes+num_dof*NUMNODESP1_2d*3+num_dof*i+1)=0; + *(LStokes+num_dof*NUMNODESP1_2d*3+num_dof*i+2)=l1l2l3[i]; + *(LStokes+num_dof*NUMNODESP1_2d*3+num_dof*i+3)=0; + + } +} +/*}}}*/ +/*FUNCTION PentaRef::GetLprimeStokesMacAyeal {{{1*/ +void PentaRef::GetLprimeStokesMacAyeal(double* LprimeStokes, double* xyz_list, GaussPenta* gauss){ + + /* + * Compute Lprime matrix. Lprime=[Lp1 Lp2 Lp3] where Lpi is square and of size numdof. + * For grid i, Lpi can be expressed in the actual coordinate system + * by: + * Lpi=[ h 0 ] + * [ 0 h ] + * [ h 0 ] + * [ 0 h ] + * where h is the interpolation function for grid i. + */ + int i; + int num_dof=2; + + double l1l2l3[NUMNODESP1_2d]; + double dh1dh6[3][NUMNODESP1]; + + /*Get l1l2l3 in actual coordinate system: */ + l1l2l3[0]=gauss->coord1*(1-gauss->coord4)/2.0; + l1l2l3[1]=gauss->coord2*(1-gauss->coord4)/2.0; + l1l2l3[2]=gauss->coord3*(1-gauss->coord4)/2.0; + + GetNodalFunctionsP1Derivatives(&dh1dh6[0][0],xyz_list,gauss); + + /*Build LprimeStokes: */ + for (i=0;i<3;i++){ + *(LprimeStokes+num_dof*NUMNODESP1_2d*0+num_dof*i)=l1l2l3[i]; //LprimeStokes[0][NDOF2*i]=dh1dh3[0][i]; + *(LprimeStokes+num_dof*NUMNODESP1_2d*0+num_dof*i+1)=0; + *(LprimeStokes+num_dof*NUMNODESP1_2d*1+num_dof*i)=0; + *(LprimeStokes+num_dof*NUMNODESP1_2d*1+num_dof*i+1)=l1l2l3[i]; + *(LprimeStokes+num_dof*NUMNODESP1_2d*2+num_dof*i)=l1l2l3[i]; + *(LprimeStokes+num_dof*NUMNODESP1_2d*2+num_dof*i+1)=0; + *(LprimeStokes+num_dof*NUMNODESP1_2d*3+num_dof*i)=0; + *(LprimeStokes+num_dof*NUMNODESP1_2d*3+num_dof*i+1)=l1l2l3[i]; + } +} +/*}}}*/ /*FUNCTION PentaRef::GetJacobian {{{1*/ void PentaRef::GetJacobian(double* J, double* xyz_list,GaussPenta* gauss){ Index: /issm/trunk/src/c/objects/Elements/PentaRef.h =================================================================== --- /issm/trunk/src/c/objects/Elements/PentaRef.h (revision 7069) +++ /issm/trunk/src/c/objects/Elements/PentaRef.h (revision 7070) @@ -51,4 +51,8 @@ void GetLStokes(double* LStokes, GaussPenta* gauss); void GetLprimeStokes(double* LprimeStokes, double* xyz_list, GaussPenta* gauss); + void GetLMacAyealStokes(double* LMacAyealStokes, GaussPenta* gauss); + void GetLprimeMacAyealStokes(double* LprimeMacAyealStokes, double* xyz_list, GaussPenta* gauss); + void GetLStokesMacAyeal(double* LStokesMacAyeal, GaussPenta* gauss); + void GetLprimeStokesMacAyeal(double* LprimeStokesMacAyeal, double* xyz_list, GaussPenta* gauss); void GetParameterValue(double* pvalue,double* plist, GaussPenta* gauss); void GetParameterValue(double* pvalue,double* plist,GaussTria* gauss){_error_("only PentaGauss are supported");};