Index: /issm/trunk/src/c/objects/Elements/Penta.cpp =================================================================== --- /issm/trunk/src/c/objects/Elements/Penta.cpp (revision 6986) +++ /issm/trunk/src/c/objects/Elements/Penta.cpp (revision 6987) @@ -703,32 +703,115 @@ /*compute all stiffness matrices for this element*/ - 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; - Ke1=CreateKMatrixDiagnosticMacAyeal3d(); - Ke2=CreateKMatrixDiagnosticStokes(); - - _error_("KMatrix coupling MacAyealStokes not terminated yet"); - - /*Constants*/ - const int numdofm=NDOF2*NUMVERTICES; - const int numdofs=NDOF4*NUMVERTICES; - const int numdoftotal=(NDOF2+NDOF4)*NUMVERTICES; - int i,j; - - for(i=0;ivalues[(i+numdofm)*numdoftotal+NDOF2*j+0]+=Ke2->values[i*numdofs+NDOF4*j+0]; - Ke->values[(i+numdofm)*numdoftotal+NDOF2*j+1]+=Ke2->values[i*numdofs+NDOF4*j+1]; - } - for(i=0;ivalues[i*numdoftotal+numdofm+NDOF4*j+0]+=Ke1->values[i*numdofm+NDOF2*j+0]; - Ke->values[i*numdoftotal+numdofm+NDOF4*j+1]+=Ke1->values[i*numdofm+NDOF2*j+1]; - } + ElementMatrix* Ke1=CreateKMatrixCouplingMacAyealPattynViscous(); + ElementMatrix* Ke2=CreateKMatrixCouplingMacAyealPattynFriction(); + ElementMatrix* Ke =new ElementMatrix(Ke1,Ke2); /*clean-up and return*/ delete Ke1; delete Ke2; + return Ke; +} +/*}}}*/ +/*FUNCTION Penta::CreateKMatrixCouplingMacAyealStokesViscous{{{1*/ +ElementMatrix* Penta::CreateKMatrixCouplingMacAyealStokesViscous(void){ + + /*Constants*/ + const int numdofm=NDOF2*NUMVERTICES2D; + const int numdofs=NDOF4*NUMVERTICES; + const int numdoftotal=numdofm+numdofs; + + /*Intermediaries */ + int i,j,ig; + double Jdet; + double viscosity; //viscosity + double epsilon[5]; /* epsilon=[exx,eyy,exy,exz,eyz];*/ + double xyz_list[NUMVERTICES][3]; + double B[4][numdofs]; + double Bprime[4][numdofm]; + double B2[3][numdofm]; + double Bprime2[3][numdofs]; + double D[4][4]={0.0}; // material matrix, simple scalar matrix. + double D2[3][3]={0.0}; // material matrix, simple scalar matrix. + double D_scalar; + double Ke_gg[numdofs][numdofm]={0.0}; //local element stiffness matrix + double Ke_gg2[numdofm][numdofs]={0.0}; //local element stiffness matrix + double Ke_gg_gaussian[numdofs][numdofm]; //stiffness matrix evaluated at the gaussian point. + double Ke_gg_gaussian2[numdofm][numdofs]; //stiffness matrix evaluated at the gaussian point. + GaussPenta *gauss=NULL; + GaussTria *gauss_tria=NULL; + + /*Find penta on bed as stokes must be coupled to the dofs on the bed: */ + Penta* pentabase=GetBasalElement(); + Tria* tria=pentabase->SpawnTria(0,1,2); //grids 0, 1 and 2 make the new tria. + + /*Initialize Element matrix and return if necessary*/ + if(IsOnWater()) return NULL; + ElementMatrix* Ke1=new ElementMatrix(pentabase->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; + + /* Get node coordinates and dof list: */ + GetVerticesCoordinates(&xyz_list[0][0], nodes, 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); + + /* Start looping on the number of gaussian points: */ + gauss=new GaussPenta(5,5); + gauss_tria=new GaussTria(); + for (ig=gauss->begin();igend();ig++){ + + gauss->GaussPoint(ig); + gauss->SynchronizeGaussTria(gauss_tria); + + GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss); + GetBMacAyealStokes(&B[0][0], &xyz_list[0][0], gauss); + tria->GetBprimeMacAyealStokes(&Bprime[0][0], &xyz_list[0][0], gauss_tria); + tria->GetBMacAyealStokes(&B2[0][0], &xyz_list[0][0], gauss_tria); + GetBprimeMacAyealStokes(&Bprime2[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]); + + D_scalar=2*viscosity*gauss->weight*Jdet; + for (i=0;i<3;i++) D[i][i]=D_scalar; + D[4][4]=(-1)*gauss->weight*Jdet; + for (i=0;i<3;i++) D2[i][i]=D_scalar; + + TripleMultiply( &B[0][0],4,numdofs,1, + &D[0][0],4,4,0, + &Bprime[0][0],4,numdofm,0, + &Ke_gg_gaussian[0][0],0); + + TripleMultiply( &B2[0][0],3,numdofm,1, + &D2[0][0],3,3,0, + &Bprime2[0][0],3,numdofs,0, + &Ke_gg_gaussian2[0][0],0); + + for( i=0; ivalues[(i+2*numdofm)*numdoftotal+j]+=Ke_gg[i][j]; + for(i=0;ivalues[i*numdoftotal+(j+2*numdofm)]+=Ke_gg2[j][i]; + _error_("coupling not finished yet"); + + /*Clean-up and return*/ + delete tria->matice; delete tria; + delete gauss; + delete gauss_tria; + return Ke; +} +/*}}}*/ +/*FUNCTION Penta::CreateKMatrixCouplingMacAyealStokesFriction{{{1*/ +ElementMatrix* Penta::CreateKMatrixCouplingMacAyealStokesFriction(void){ + + /*Initialize Element matrix and return if necessary*/ + 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; + return Ke; } @@ -903,5 +986,5 @@ /*Intermediaries */ - int i,j,ig; + int i,j,ig,approximation; double Jdet; double viscosity, oldviscosity, newviscosity, viscosity_overshoot; @@ -925,4 +1008,5 @@ if(IsOnWater()) return NULL; ElementMatrix* Ke=new ElementMatrix(tria->nodes,NUMVERTICES2D,this->parameters,MacAyealApproximationEnum); + inputs->GetParameterValue(&approximation,ApproximationEnum); /*Retrieve all inputs and parameters*/ @@ -933,4 +1017,5 @@ Input* vxold_input=inputs->GetInput(VxOldEnum); _assert_(vxold_input); Input* vyold_input=inputs->GetInput(VyOldEnum); _assert_(vyold_input); + Input* vz_input=inputs->GetInput(VzEnum); _assert_(vz_input); /* Start looping on the number of gaussian points: */ @@ -946,10 +1031,18 @@ tria->GetBprimeMacAyeal(&Bprime[0][0], &xyz_list[0][0], gauss_tria); - this->GetStrainRate3dPattyn(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input); - this->GetStrainRate3dPattyn(&oldepsilon[0],&xyz_list[0][0],gauss,vxold_input,vyold_input); - matice->GetViscosity3d(&viscosity, &epsilon[0]); - matice->GetViscosity3d(&oldviscosity, &oldepsilon[0]); - - newviscosity=viscosity+viscosity_overshoot*(viscosity-oldviscosity); + if(approximation==MacAyealPattynApproximationEnum){ + this->GetStrainRate3dPattyn(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input); + this->GetStrainRate3dPattyn(&oldepsilon[0],&xyz_list[0][0],gauss,vxold_input,vyold_input); + matice->GetViscosity3d(&viscosity, &epsilon[0]); + matice->GetViscosity3d(&oldviscosity, &oldepsilon[0]); + + newviscosity=viscosity+viscosity_overshoot*(viscosity-oldviscosity); + } + else if (approximation==MacAyealStokesApproximationEnum){ + this->GetStrainRate3d(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input,vz_input); + matice->GetViscosity3dStokes(&newviscosity,&epsilon[0]); + } + else _error_("approximation %i (%s) not supported yet",approximation,EnumToString(approximation)); + D_scalar=2*newviscosity*gauss->weight*Jdet; for (i=0;i<3;i++) D[i][i]=D_scalar; Index: /issm/trunk/src/c/objects/Elements/Penta.h =================================================================== --- /issm/trunk/src/c/objects/Elements/Penta.h (revision 6986) +++ /issm/trunk/src/c/objects/Elements/Penta.h (revision 6987) @@ -131,4 +131,6 @@ ElementMatrix* CreateKMatrixCouplingMacAyealPattynFriction(void); ElementMatrix* CreateKMatrixCouplingMacAyealStokes(void); + ElementMatrix* CreateKMatrixCouplingMacAyealStokesViscous(void); + ElementMatrix* CreateKMatrixCouplingMacAyealStokesFriction(void); ElementMatrix* CreateKMatrixCouplingPattynStokes(void); ElementMatrix* CreateKMatrixDiagnosticHoriz(void); Index: /issm/trunk/src/c/objects/Elements/PentaRef.cpp =================================================================== --- /issm/trunk/src/c/objects/Elements/PentaRef.cpp (revision 6986) +++ /issm/trunk/src/c/objects/Elements/PentaRef.cpp (revision 6987) @@ -86,4 +86,51 @@ *(B+NDOF2*NUMNODESP1*2+NDOF2*i+1)=(float).5*dh1dh6[0][i]; + } + +} +/*}}}*/ +/*FUNCTION PentaRef::GetBMacAyealStokes{{{1*/ +void PentaRef::GetBMacAyealStokes(double* B, double* xyz_list, GaussPenta* gauss){ + /*Compute B matrix. B=[B1 B2 B3 B4 B5 B6] where Bi is of size 5*NDOF2. + * For grid i, Bi can be expressed in the actual coordinate system + * by: + * Bi=[ dh/dx 0 0 0 ] + * [ 0 dh/dy 0 0 ] + * [ 1/2*dh/dy 1/2*dh/dx 0 0 ] + * [ 0 0 0 h ] + * where h is the interpolation function for grid i. + * + * We assume B has been allocated already, of size: 5x(NDOF2*NUMNODESP1) + */ + + int i; + double dh1dh7[3][NUMNODESMINI]; + double l1l6[NUMNODESP1]; + + /*Get dh1dh6 in actual coordinate system: */ + GetNodalFunctionsMINIDerivatives(&dh1dh7[0][0],xyz_list, gauss); + GetNodalFunctionsP1(l1l6, gauss); + + /*Build B: */ + for (i=0;i