Index: /issm/trunk/src/c/objects/Elements/Penta.cpp =================================================================== --- /issm/trunk/src/c/objects/Elements/Penta.cpp (revision 5653) +++ /issm/trunk/src/c/objects/Elements/Penta.cpp (revision 5654) @@ -2607,43 +2607,185 @@ void Penta::CreateKMatrixDiagnosticMacAyeal( Mat Kgg){ + /* local declarations */ + int i,j; + + /* node data: */ + const int numgrids3d=6; + const int numgrids2d=3; + const int numdof2d=2*numgrids2d; + double xyz_list[numgrids3d][3]; + int* doflist=NULL; + + /* 3d gaussian points: */ + int ig; + GaussPenta *gauss=NULL; + GaussTria *gauss_tria=NULL; + + /* material data: */ + double viscosity; //viscosity + double oldviscosity; //viscosity + double newviscosity; //viscosity + + /* strain rate: */ + double epsilon[5]; /* epsilon=[exx,eyy,exy,exz,eyz];*/ + double oldepsilon[5]; /* epsilon=[exx,eyy,exy,exz,eyz];*/ + + /* matrices: */ + double B[3][numdof2d]; + double Bprime[3][numdof2d]; + double L[2][numdof2d]; + double D[3][3]={0.0}; // material matrix, simple scalar matrix. + double D_scalar; + double DL[2][2]={0.0}; //for basal drag + double DL_scalar; + + /* local element matrices: */ + double Ke_gg[numdof2d][numdof2d]={0.0}; //local element stiffness matrix + double Ke_gg_gaussian[numdof2d][numdof2d]; //stiffness matrix evaluated at the gaussian point. + double Jdet; + + /*slope: */ + double slope[2]={0.0}; + double slope_magnitude; + + /*friction: */ + double alpha2_list[3]; + double alpha2; + + double MAXSLOPE=.06; // 6 % + double MOUNTAINKEXPONENT=10; + + /*parameters: */ + double viscosity_overshoot; + /*Collapsed formulation: */ Tria* tria=NULL; + Penta* pentabase=NULL; /*inputs: */ bool onwater; bool onbed; - + bool shelf; + int approximation; + Input* vx_input=NULL; + Input* vy_input=NULL; + Input* vxold_input=NULL; + Input* vyold_input=NULL; + + inputs->GetParameterValue(&approximation,ApproximationEnum); inputs->GetParameterValue(&onwater,ElementOnWaterEnum); inputs->GetParameterValue(&onbed,ElementOnBedEnum); + inputs->GetParameterValue(&shelf,ElementOnIceShelfEnum); /*If on water, skip stiffness: */ - if(onwater)return; - - /*Figure out if this pentaelem is collapsed. If so, then bailout, except if it is at the - bedrock, in which case we spawn a tria element using the 3 first grids, and use it to build - the stiffness matrix. */ - - if (onbed==0){ - /*This element should be collapsed, but this element is not on the bedrock, therefore all its - * dofs have already been frozen! Do nothing: */ - return; - } - else if (onbed==1){ - - /*This element should be collapsed into a tria element at its base. Create this tria element, - *and use its CreateKMatrix functionality to fill the global stiffness matrix: */ - - /*Depth Averaging B*/ - this->InputDepthAverageAtBase(RheologyBEnum,RheologyBbarEnum,MaterialsEnum); - - /*Call Tria function*/ - tria=(Tria*)SpawnTria(0,1,2); //grids 0, 1 and 2 make the new tria. - tria->CreateKMatrix(Kgg); - delete tria->matice; delete tria; - - /*Delete B averaged*/ - this->matice->inputs->DeleteInput(RheologyBbarEnum); - - return; + if(approximation==MacAyealApproximationEnum){ + if(onwater)return; + + /*Figure out if this pentaelem is collapsed. If so, then bailout, except if it is at the + bedrock, in which case we spawn a tria element using the 3 first grids, and use it to build + the stiffness matrix. */ + + if (onbed==0){ + /*This element should be collapsed, but this element is not on the bedrock, therefore all its + * dofs have already been frozen! Do nothing: */ + return; + } + else if (onbed==1){ + + /*This element should be collapsed into a tria element at its base. Create this tria element, + *and use its CreateKMatrix functionality to fill the global stiffness matrix: */ + + /*Depth Averaging B*/ + this->InputDepthAverageAtBase(RheologyBEnum,RheologyBbarEnum,MaterialsEnum); + + /*Call Tria function*/ + tria=(Tria*)SpawnTria(0,1,2); //grids 0, 1 and 2 make the new tria. + tria->CreateKMatrix(Kgg); + delete tria->matice; delete tria; + + /*Delete B averaged*/ + this->matice->inputs->DeleteInput(RheologyBbarEnum); + + return; + } + } + else if(approximation==MacAyealPattynApproximationEnum){ + /*retrieve some parameters: */ + this->parameters->FindParam(&viscosity_overshoot,ViscosityOvershootEnum); + + /*If on water, skip stiffness: */ + if(onwater)return; + + /*Find penta on bed as this is a macayeal elements: */ + pentabase=GetBasalElement(); + tria=pentabase->SpawnTria(0,1,2); //grids 0, 1 and 2 make the new tria. + + /* Get node coordinates and dof list: */ + GetVerticesCoordinates(&xyz_list[0][0], nodes, numgrids3d); + tria->GetDofList(&doflist,MacAyealApproximationEnum); //Pattyn dof list + + /*Retrieve all inputs we will be needing: */ + vx_input=inputs->GetInput(VxEnum); + vy_input=inputs->GetInput(VyEnum); + vxold_input=inputs->GetInput(VxOldEnum); + vyold_input=inputs->GetInput(VyOldEnum); + + /* 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); + + /*Get strain rate from velocity: */ + 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); + + /*Get viscosity: */ + matice->GetViscosity3d(&viscosity, &epsilon[0]); + matice->GetViscosity3d(&oldviscosity, &oldepsilon[0]); + + /*Get B and Bprime matrices: */ + tria->GetBMacAyeal(&B[0][0], &xyz_list[0][0], gauss_tria); + tria->GetBprimeMacAyeal(&Bprime[0][0], &xyz_list[0][0], gauss_tria); + + /* Get Jacobian determinant: */ + GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss); + + /*Build the D matrix: we plug the gaussian weight, the viscosity, and the jacobian determinant + onto this scalar matrix, so that we win some computational time: */ + newviscosity=viscosity+viscosity_overshoot*(viscosity-oldviscosity); + D_scalar=2*newviscosity*gauss->weight*Jdet; + for (i=0;i<3;i++) D[i][i]=D_scalar; + + /* Do the triple product tB*D*Bprime: */ + TripleMultiply( &B[0][0],3,numdof2d,1, + &D[0][0],3,3,0, + &Bprime[0][0],3,numdof2d,0, + &Ke_gg_gaussian[0][0],0); + + /* Add the Ke_gg_gaussian, and optionally Ke_gg_gaussian onto Ke_gg: */ + for(i=0;iCreateKMatrixDiagnosticMacAyealFriction(Kgg); + } + + /*Clean up and return*/ + delete tria->matice; + delete tria; + delete gauss_tria; + delete gauss; + xfree((void**)&doflist); } }