Index: /issm/trunk/src/c/objects/Elements/Penta.cpp =================================================================== --- /issm/trunk/src/c/objects/Elements/Penta.cpp (revision 5202) +++ /issm/trunk/src/c/objects/Elements/Penta.cpp (revision 5203) @@ -755,4 +755,10 @@ return; } + else if(approximation==MacAyealPattynApproximationEnum){ + ISSMERROR("coupling MacAyeal/Pattyn not implemented yet"); + CreateKMatrixDiagnosticMacAyeal( Kgg); + CreateKMatrixDiagnosticPattyn( Kgg); + CreateKMatrixCouplingMacAyealPattyn( Kgg); + } else ISSMERROR("Approximation %s not supported yet",EnumToString(approximation)); } @@ -827,4 +833,8 @@ else if(approximation==HutterApproximationEnum){ return; + } + else if(approximation==MacAyealPattynApproximationEnum){ + CreatePVectorDiagnosticMacAyeal( pg); + CreatePVectorDiagnosticPattyn( pg); } else ISSMERROR("Approximation %s not supported yet",EnumToString(approximation)); @@ -2083,4 +2093,207 @@ } /*}}}*/ +/*FUNCTION Penta::CreateKMatrixCouplingMacAyealPattyn{{{1*/ +void Penta::CreateKMatrixCouplingMacAyealPattyn( Mat Kgg){ + + /* local declarations */ + int i,j; + + /* node data: */ + const int numgrids=6; // Pattyn numgrids + const int numdof=2*numgrids; + const int numgrids2=3; //MacAyeal numgrids + const int numdof2=2*numgrids2; + double xyz_list[numgrids][3]; + int* doflist=NULL; + int* doflist2=NULL; + + /* 3d gaussian points: */ + int num_gauss,ig; + double* first_gauss_area_coord = NULL; + double* second_gauss_area_coord = NULL; + double* third_gauss_area_coord = NULL; + double* fourth_gauss_vert_coord = NULL; + double* area_gauss_weights = NULL; + double* vert_gauss_weights = NULL; + int ig1,ig2; + double gauss_weight1,gauss_weight2; + double gauss_coord[4]; + int order_area_gauss; + int num_vert_gauss; + int num_area_gauss; + double gauss_weight; + + /* 2d gaussian point: */ + int num_gauss2d; + double* first_gauss_area_coord2d = NULL; + double* second_gauss_area_coord2d = NULL; + double* third_gauss_area_coord2d = NULL; + double* gauss_weights2d=NULL; + double gauss_l1l2l3[3]; + + /* 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][numdof]; + double Bprime[3][numdof2]; + double L[2][numdof]; + 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[numdof][numdof2]={0.0}; //local element stiffness matrix + double Ke_gg_gaussian[numdof][numdof2]; //stiffness matrix evaluated at the gaussian point. + double Jdet; + + /*friction: */ + double alpha2_list[3]; + double alpha2; + + /*parameters: */ + double viscosity_overshoot; + + /*Collapsed formulation: */ + Tria* tria =NULL; + Penta* pentabase=NULL; + + /*inputs: */ + bool onwater; + bool onbed; + bool shelf; + Input* vx_input=NULL; + Input* vy_input=NULL; + Input* vxold_input=NULL; + Input* vyold_input=NULL; + + /*retrieve inputs :*/ + inputs->GetParameterValue(&onwater,ElementOnWaterEnum); + inputs->GetParameterValue(&onbed,ElementOnBedEnum); + inputs->GetParameterValue(&shelf,ElementOnIceShelfEnum); + + /*retrieve some parameters: */ + this->parameters->FindParam(&viscosity_overshoot,ViscosityOvershootEnum); + + /*If on water, skip stiffness: */ + if(onwater)return; + + /*Find penta on bed as pattyn must be coupled to the dofs on the bed: */ + 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, numgrids); + GetDofList(&doflist); //Pattyn dof list + GetDofList(&doflist2); //MacAyeal dof list + + /*Get gaussian points and weights. Penta is an extrusion of a Tria, we therefore + get tria gaussian points as well as segment gaussian points. For tria gaussian + points, order of integration is 2, because we need to integrate the product tB*D*B' + which is a polynomial of degree 3 (see GaussTria for more details). For segment gaussian + points, same deal, which yields 3 gaussian points.*/ + + order_area_gauss=5; + num_vert_gauss=5; + + GaussPenta( &num_area_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &area_gauss_weights, &fourth_gauss_vert_coord,&vert_gauss_weights,order_area_gauss,num_vert_gauss); + + /*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: */ + for (ig1=0; ig1GetStrainRate3dPattyn(&epsilon[0],&xyz_list[0][0],gauss_coord,vx_input,vy_input); + this->GetStrainRate3dPattyn(&oldepsilon[0],&xyz_list[0][0],gauss_coord,vxold_input,vyold_input); + + /*Get viscosity: */ + matice->GetViscosity3d(&viscosity, &epsilon[0]); + matice->GetViscosity3d(&oldviscosity, &oldepsilon[0]); + + /*Get B and Bprime matrices: */ + GetBMacAyealPattyn(&B[0][0], &xyz_list[0][0], gauss_coord); + tria->GetBprimeMacAyeal(&Bprime[0][0], &xyz_list[0][0], gauss_coord); + + /* Get Jacobian determinant: */ + GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss_coord); + + /*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=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,numdof,1, + &D[0][0],3,3,0, + &Bprime[0][0],3,numdof2,0, + &Ke_gg_gaussian[0][0],0); + + /* Add the Ke_gg_gaussian, and optionally Ke_gg_gaussian onto Ke_gg: */ + for( i=0; iCreateKMatrixDiagnosticHorizFriction(Kgg); + } + + delete tria; + delete pentabase; + + xfree((void**)&first_gauss_area_coord); + xfree((void**)&second_gauss_area_coord); + xfree((void**)&third_gauss_area_coord); + xfree((void**)&fourth_gauss_vert_coord); + xfree((void**)&area_gauss_weights); + xfree((void**)&vert_gauss_weights); + xfree((void**)&first_gauss_area_coord2d); + xfree((void**)&second_gauss_area_coord2d); + xfree((void**)&third_gauss_area_coord2d); + xfree((void**)&gauss_weights2d); + xfree((void**)&doflist); +} +/*}}}*/ /*FUNCTION Penta::CreateKMatrixDiagnosticHutter{{{1*/ void Penta::CreateKMatrixDiagnosticHutter(Mat Kgg){ @@ -2268,5 +2481,4 @@ double Ke_gg[numdof][numdof]={0.0}; //local element stiffness matrix double Ke_gg_gaussian[numdof][numdof]; //stiffness matrix evaluated at the gaussian point. - double Ke_gg_drag_gaussian[numdof][numdof]; //stiffness matrix contribution from drag double Jdet; @@ -2798,5 +3010,5 @@ &Ke_gg_gaussian[0][0],0); - /* Add the Ke_gg_gaussian, and optionally Ke_gg_drag_gaussian onto Ke_gg: */ + /* Add the Ke_gg_gaussian, and optionally Ke_gg_gaussian onto Ke_gg: */ for( i=0; i