Index: /issm/trunk/src/c/objects/Elements/Penta.cpp =================================================================== --- /issm/trunk/src/c/objects/Elements/Penta.cpp (revision 10479) +++ /issm/trunk/src/c/objects/Elements/Penta.cpp (revision 10480) @@ -5986,4 +5986,5 @@ ElementMatrix* Ke1=CreateKMatrixDiagnosticStokesViscous(); ElementMatrix* Ke2=CreateKMatrixDiagnosticStokesFriction(); + //ElementMatrix* Ke2=CreateKMatrixDiagnosticStokesFriction2(); ElementMatrix* Ke =new ElementMatrix(Ke1,Ke2); @@ -6136,9 +6137,111 @@ /*Transform Coordinate System*/ + if(id==1) printarray(&bed_normal[0],1,3); + if(id==1) printarray(&nodes[0]->coord_system[0][0],3,3); + if(id==1) printf("Before rotation\n"); + if(id==1) printarray(Ke->values,4,Ke->ncols); TransformStiffnessMatrixCoord(Ke,nodes,NUMVERTICES,NDOF4); + if(id==1) printf("After rotation\n"); + if(id==1) printarray(Ke->values,4,Ke->ncols); + if(id==1) printf("New:\n"); + if(id==1) Ke=CreateKMatrixDiagnosticStokesFriction2(); + + //if(id==56) printarray(Ke->values,1,9); + //if(id==56) printarray(&nodes[0]->coord_system[0][0],3,3); + //if(id==56) _error_("STOP"); /*Clean up and return*/ delete gauss; delete friction; + return Ke; +} +/*}}}*/ +/*FUNCTION Penta::CreateKMatrixDiagnosticStokesFriction2{{{1*/ +ElementMatrix* Penta::CreateKMatrixDiagnosticStokesFriction2(void){ + + /*Constants*/ + const int numdof=NUMVERTICES*NDOF4; + const int numdof2d=NUMVERTICES2D*NDOF4; + + /*Intermediaries */ + int i,j,ig; + int analysis_type,approximation; + double alpha2,Jdet2d; + double stokesreconditioning,viscosity; + double epsilon[6]; /* epsilon=[exx,eyy,ezz,exy,exz,eyz];*/ + double xyz_list[NUMVERTICES][3]; + double xyz_list_tria[NUMVERTICES2D][3]; + double LStokes2[2][numdof2d]; + //double LprimeStokes2[4][numdof2d]; + double DLStokes[2][2]={0.0}; + double 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], nodes, 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); + GetLStokes3(&LStokes2[0][0], gauss); + //GetLprimeStokes2(&LprimeStokes2[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]); + + 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 + //DLStokes[2][2] = +2*viscosity*gauss->weight*Jdet2d; + //DLStokes[3][3] = -stokesreconditioning*gauss->weight*Jdet2d; + + //TripleMultiply( &LStokes2[0][0],4,numdof2d,1, + // &DLStokes[0][0],4,4,0, + // &LprimeStokes2[0][0],4,numdof2d,0, + // &Ke_drag_gaussian[0][0],0); + TripleMultiply( &LStokes2[0][0],2,numdof2d,1, + &DLStokes[0][0],2,2,0, + &LStokes2[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,NDOF4); + + //printf("New:\n"); + //printarray(&Ke_drag_gaussian[0][0],numdof2d,numdof2d); + //printf("Old:\n"); + //Ke=CreateKMatrixDiagnosticStokesFriction(); + //_error_("STOP"); + //if(id==1) printarray(Ke->values,4,Ke->ncols); + //if(id==1) printarray(Ke->values,1,9); + + /*Clean up and return*/ + delete gauss; + delete friction; + + //return CreateKMatrixDiagnosticStokesFriction(); return Ke; } Index: /issm/trunk/src/c/objects/Elements/Penta.h =================================================================== --- /issm/trunk/src/c/objects/Elements/Penta.h (revision 10479) +++ /issm/trunk/src/c/objects/Elements/Penta.h (revision 10480) @@ -223,4 +223,5 @@ ElementMatrix* CreateKMatrixDiagnosticStokesViscous(void); ElementMatrix* CreateKMatrixDiagnosticStokesFriction(void); + ElementMatrix* CreateKMatrixDiagnosticStokesFriction2(void); ElementMatrix* CreateKMatrixDiagnosticVert(void); ElementMatrix* CreateKMatrixDiagnosticVertVolume(void); Index: /issm/trunk/src/c/objects/Elements/PentaRef.cpp =================================================================== --- /issm/trunk/src/c/objects/Elements/PentaRef.cpp (revision 10479) +++ /issm/trunk/src/c/objects/Elements/PentaRef.cpp (revision 10480) @@ -651,4 +651,141 @@ } /*}}}*/ +/*FUNCTION PentaRef::GetLStokes2{{{1*/ +void PentaRef::GetLStokes2(double* LStokes, GaussPenta* gauss){ + /* + * Compute L matrix. L=[L1 L2 L3] where Li is square and of size numdof. + * For node i, Li can be expressed in the actual coordinate system + * by: + * Li=[ h 0 0 0] + * [ 0 h 0 0] + * [ 0 0 h h] + * [ 0 0 0 0] + * where h is the interpolation function for node 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+0)=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)=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)=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)=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::GetLStokes3{{{1*/ +void PentaRef::GetLStokes3(double* LStokes, GaussPenta* gauss){ + /* + * Compute L matrix. L=[L1 L2 L3] where Li is square and of size numdof. + * For node i, Li can be expressed in the actual coordinate system + * by: + * Li=[ h 0 ] + * [ 0 h ] + * [ 0 0 ] + * [ 0 0 ] + * where h is the interpolation function for node 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+0)=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)=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.; + } +} +/*}}}*/ +/*FUNCTION PentaRef::GetLprimeStokes2{{{1*/ +void PentaRef::GetLprimeStokes2(double* LStokes,double* xyz_list,GaussPenta* gauss){ + /* + * Compute L matrix. L=[L1 L2 L3] where Li is square and of size numdof. + * For node i, Li can be expressed in the actual coordinate system + * by: + * Li=[ h 0 0 0] + * [ 0 h 0 0] + * [ 0 0 dh\dz 0] + * [ 0 0 0 h] + * where h is the interpolation function for node 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 LStokes: */ + for (i=0;i<3;i++){ + *(LStokes+num_dof*NUMNODESP1_2d*0+num_dof*i+0)=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)=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)=0.; + *(LStokes+num_dof*NUMNODESP1_2d*2+num_dof*i+1)=0.; + *(LStokes+num_dof*NUMNODESP1_2d*2+num_dof*i+2)=dh1dh6[2][i]; + *(LStokes+num_dof*NUMNODESP1_2d*2+num_dof*i+3)=0.; + + *(LStokes+num_dof*NUMNODESP1_2d*3+num_dof*i+0)=0.; + *(LStokes+num_dof*NUMNODESP1_2d*3+num_dof*i+1)=0.; + *(LStokes+num_dof*NUMNODESP1_2d*3+num_dof*i+2)=0.; + *(LStokes+num_dof*NUMNODESP1_2d*3+num_dof*i+3)=l1l2l3[i]; + } +} +/*}}}*/ /*FUNCTION PentaRef::GetLprimeStokes {{{1*/ void PentaRef::GetLprimeStokes(double* LprimeStokes, double* xyz_list, GaussPenta* gauss){ @@ -658,18 +795,33 @@ * For node i, Lpi can be expressed in the actual coordinate system * by: - * Lpi=[ h 0 0 0] - * [ 0 h 0 0] - * [ 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 dh/dz 0] - * [dh/dz 0 dh/dx 0] - * [ 0 dh/dz dh/dy 0] - * [ 0 0 0 h] - * [ 0 0 0 h] - * [ 0 0 0 h] + * Lpi=[ h 0 0 0]1 + * [ 0 h 0 0]2 + * [ h 0 0 0]3 + * [ 0 h 0 0]4 + * [ 0 0 h 0]5 + * [ 0 0 h 0]6 + * [ 0 0 dh/dz 0]7 + * [ 0 0 dh/dz 0]8 + * [ 0 0 dh/dz 0]9 + * [dh/dz 0 dh/dx 0]0 + * [ 0 dh/dz dh/dy 0]1 + * [ 0 0 0 h]2 + * [ 0 0 0 h]3 + * [ 0 0 0 h]4 + * + * Li=[ h 0 0 0]1 + * [ 0 h 0 0]2 + * [ 0 0 h 0]3 + * [ 0 0 h 0]4 + * [ h 0 0 0]5 + * [ 0 h 0 0]6 + * [ h 0 0 0]7 + * [ 0 h 0 0]8 + * [ 0 0 h 0]9 + * [ 0 0 h 0]0 + * [ 0 0 h 0]1 + * [ h 0 0 0]2 + * [ 0 h 0 0]3 + * [ 0 0 h 0]4 * where h is the interpolation function for node i. */ Index: /issm/trunk/src/c/objects/Elements/PentaRef.h =================================================================== --- /issm/trunk/src/c/objects/Elements/PentaRef.h (revision 10479) +++ /issm/trunk/src/c/objects/Elements/PentaRef.h (revision 10480) @@ -51,4 +51,7 @@ void GetL(double* L, GaussPenta* gauss,int numdof); void GetLStokes(double* LStokes, GaussPenta* gauss); + void GetLStokes2(double* LStokes, GaussPenta* gauss); + void GetLStokes3(double* LStokes, GaussPenta* gauss); + void GetLprimeStokes2(double* LStokes,double* xyz_list,GaussPenta* gauss); void GetLprimeStokes(double* LprimeStokes, double* xyz_list, GaussPenta* gauss); void GetLMacAyealStokes(double* LMacAyealStokes, GaussPenta* gauss); Index: /issm/trunk/src/c/shared/Elements/elements.h =================================================================== --- /issm/trunk/src/c/shared/Elements/elements.h (revision 10479) +++ /issm/trunk/src/c/shared/Elements/elements.h (revision 10480) @@ -28,5 +28,5 @@ for(int i=0;i