Index: /issm/trunk/src/c/objects/Elements/Penta.cpp =================================================================== --- /issm/trunk/src/c/objects/Elements/Penta.cpp (revision 5846) +++ /issm/trunk/src/c/objects/Elements/Penta.cpp (revision 5847) @@ -2193,5 +2193,5 @@ switch(approximation){ case MacAyealApproximationEnum: - CreateKMatrixDiagnosticMacAyeal( Kgg); + CreateKMatrixDiagnosticMacAyeal2d( Kgg); break; case PattynApproximationEnum: @@ -2206,5 +2206,5 @@ break; case MacAyealPattynApproximationEnum: - CreateKMatrixDiagnosticMacAyeal( Kgg); + CreateKMatrixDiagnosticMacAyeal3d( Kgg); CreateKMatrixDiagnosticPattyn( Kgg); CreateKMatrixCouplingMacAyealPattyn( Kgg); @@ -2290,30 +2290,51 @@ } -/*FUNCTION Penta::CreateKMatrixDiagnosticMacAyeal{{{1*/ -void Penta::CreateKMatrixDiagnosticMacAyeal( Mat Kgg){ +/*FUNCTION Penta::CreateKMatrixDiagnosticMacAyeal2d{{{1*/ +void Penta::CreateKMatrixDiagnosticMacAyeal2d(Mat Kgg){ - /* local declarations */ - int i,j; - - /* node data: */ + /*If on water, skip stiffness: */ + if(IsOnWater()) 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 (!IsOnBed()){ + /*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{ + + /*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=(Tria*)SpawnTria(0,1,2); //grids 0, 1 and 2 make the new tria. + tria->CreateKMatrix(Kgg,NULL,NULL); + delete tria->matice; delete tria; + + /*Delete B averaged*/ + this->matice->inputs->DeleteInput(RheologyBbarEnum); + + return; + } +} +/*}}}*/ +/*FUNCTION Penta::CreateKMatrixDiagnosticMacAyeal3d{{{1*/ +void Penta::CreateKMatrixDiagnosticMacAyeal3d( Mat Kgg){ + const int numdof2d=2*NUMVERTICES2D; + int i,j,ig; double xyz_list[NUMVERTICES][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 viscosity, oldviscosity, newviscosity, viscosity_overshoot; 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]; @@ -2323,142 +2344,80 @@ 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; - int approximation; - - inputs->GetParameterValue(&approximation,ApproximationEnum); /*If on water, skip stiffness: */ - if(approximation==MacAyealApproximationEnum){ - if(IsOnWater())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 (!IsOnBed()){ - /*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{ - - /*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,NULL,NULL); - 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(IsOnWater())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,NUMVERTICES); - tria->GetDofList(&doflist,MacAyealApproximationEnum,GsetEnum); //Pattyn dof list - - /*Retrieve all inputs we will be needing: */ - Input* vx_input=inputs->GetInput(VxEnum); ISSMASSERT(vx_input); - Input* vy_input=inputs->GetInput(VyEnum); ISSMASSERT(vy_input); - Input* vxold_input=inputs->GetInput(VxOldEnum); ISSMASSERT(vxold_input); - Input* vyold_input=inputs->GetInput(VyOldEnum); ISSMASSERT(vyold_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); - - /*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,NULL,NULL); - } - - /*Clean up and return*/ - delete tria->matice; - delete tria; - delete gauss_tria; - delete gauss; - xfree((void**)&doflist); - } + if(IsOnWater())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 dof list: */ + tria->GetDofList(&doflist,MacAyealApproximationEnum,GsetEnum); //Pattyn dof list + + /*Retrieve all inputs and parameters*/ + GetVerticesCoordinates(&xyz_list[0][0], nodes,NUMVERTICES); + this->parameters->FindParam(&viscosity_overshoot,ViscosityOvershootEnum); + Input* vx_input=inputs->GetInput(VxEnum); ISSMASSERT(vx_input); + Input* vy_input=inputs->GetInput(VyEnum); ISSMASSERT(vy_input); + Input* vxold_input=inputs->GetInput(VxOldEnum); ISSMASSERT(vxold_input); + Input* vyold_input=inputs->GetInput(VyOldEnum); ISSMASSERT(vyold_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); + + tria->GetBMacAyeal(&B[0][0], &xyz_list[0][0], gauss_tria); + 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); + GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss); + D_scalar=2*newviscosity*gauss->weight*Jdet; + for (i=0;i<3;i++) D[i][i]=D_scalar; + + 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); + + for(i=0;iCreateKMatrixDiagnosticMacAyealFriction(Kgg,NULL,NULL); + } + + /*Clean up and return*/ + delete tria->matice; + delete tria; + delete gauss_tria; + delete gauss; + xfree((void**)&doflist); } /*}}}*/ Index: /issm/trunk/src/c/objects/Elements/Penta.h =================================================================== --- /issm/trunk/src/c/objects/Elements/Penta.h (revision 5846) +++ /issm/trunk/src/c/objects/Elements/Penta.h (revision 5847) @@ -127,5 +127,6 @@ void CreateKMatrixDiagnosticHoriz( Mat Kgg); void CreateKMatrixDiagnosticHutter( Mat Kgg); - void CreateKMatrixDiagnosticMacAyeal( Mat Kgg); + void CreateKMatrixDiagnosticMacAyeal2d(Mat Kgg); + void CreateKMatrixDiagnosticMacAyeal3d(Mat Kgg); void CreateKMatrixDiagnosticPattyn( Mat Kgg); void CreateKMatrixDiagnosticStokes( Mat Kgg);