Index: /issm/trunk/src/c/objects/Elements/Penta.cpp =================================================================== --- /issm/trunk/src/c/objects/Elements/Penta.cpp (revision 5871) +++ /issm/trunk/src/c/objects/Elements/Penta.cpp (revision 5872) @@ -2004,40 +2004,34 @@ /*}}}*/ /*FUNCTION Penta::CreateKMatrixCouplingMacAyealPattyn{{{1*/ -void Penta::CreateKMatrixCouplingMacAyealPattyn( Mat Kgg){ - - this->CreateKMatrixCouplingMacAyealPattynViscous(Kgg); - this->CreateKMatrixCouplingMacAyealPattynFriction(Kgg); +ElementMatrix* Penta::CreateKMatrixCouplingMacAyealPattyn(void){ + + /*compute all stiffness matrices for this element*/ + ElementMatrix* Ke1=CreateKMatrixCouplingMacAyealPattynViscous(); + ElementMatrix* Ke2=CreateKMatrixCouplingMacAyealPattynFriction(); + ElementMatrix* Ke =new ElementMatrix(Ke1,Ke2); + + /*clean-up and return*/ + delete Ke1; + delete Ke2; + return Ke; } /*}}}*/ /*FUNCTION Penta::CreateKMatrixCouplingMacAyealPattynViscous{{{1*/ -void Penta::CreateKMatrixCouplingMacAyealPattynViscous( Mat Kgg){ - - /* local declarations */ - int i,j; - - /* node data: */ +ElementMatrix* Penta::CreateKMatrixCouplingMacAyealPattynViscous(void){ + const int NUMVERTICESm=3; //MacAyealNUMVERTICES const int numdofm=2*NUMVERTICESm; const int NUMVERTICESp=6; //Pattyn NUMVERTICES const int numdofp=2*NUMVERTICESp; + const int numdoftotal=2*NDOF2*NUMVERTICES; + int i,j,ig; double xyz_list[NUMVERTICESp][3]; - int* doflistm=NULL; - int* doflistp=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][numdofp]; double Bprime[3][numdofm]; @@ -2047,38 +2041,25 @@ double DL[2][2]={0.0}; //for basal drag double DL_scalar; - - /* local element matrices: */ double Ke_gg[numdofp][numdofm]={0.0}; //local element stiffness matrix - double Ke_gg_transp[numdofm][numdofp]={0.0}; //local element stiffness matrix double Ke_gg_gaussian[numdofp][numdofm]; //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; - - /*retrieve some parameters: */ - this->parameters->FindParam(&viscosity_overshoot,ViscosityOvershootEnum); - - /*If on water, skip stiffness: */ - if(IsOnWater())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. + 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=pentabase->NewElementMatrix(MacAyealApproximationEnum); + ElementMatrix* Ke2=this->NewElementMatrix(PattynApproximationEnum); + ElementMatrix* Ke=new ElementMatrix(Ke1,Ke2); + delete Ke1; delete Ke2; /* Get node coordinates and dof list: */ GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICESp); - tria->GetDofList(&doflistm,MacAyealApproximationEnum,GsetEnum); //Pattyn dof list - GetDofList(&doflistp,PattynApproximationEnum,GsetEnum); //MacAyeal dof list - - /*Retrieve all inputs we will be needing: */ + this->parameters->FindParam(&viscosity_overshoot,ViscosityOvershootEnum); Input* vx_input=inputs->GetInput(VxEnum); ISSMASSERT(vx_input); Input* vy_input=inputs->GetInput(VyEnum); ISSMASSERT(vy_input); @@ -2094,27 +2075,17 @@ 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: */ GetBMacAyealPattyn(&B[0][0], &xyz_list[0][0], gauss); 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,numdofp,1, &D[0][0],3,3,0, @@ -2122,35 +2093,31 @@ &Ke_gg_gaussian[0][0],0); - /* Add the Ke_gg_gaussian, and optionally Ke_gg_gaussian onto Ke_gg: */ for( i=0; ivalues[(i+2*numdofm)*numdoftotal+j]+=Ke_gg[i][j]; + for(i=0;ivalues[i*numdoftotal+(j+2*numdofm)]+=Ke_gg[j][i]; + + /*Clean-up and return*/ delete tria->matice; delete tria; delete gauss; delete gauss_tria; + return Ke; } /*}}}*/ /*FUNCTION Penta::CreateKMatrixCouplingMacAyealPattynFriction{{{1*/ -void Penta::CreateKMatrixCouplingMacAyealPattynFriction( Mat Kgg){ +ElementMatrix* Penta::CreateKMatrixCouplingMacAyealPattynFriction(void){ /*Initialize Element matrix and return if necessary*/ - if(IsOnWater() || IsOnShelf() || !IsOnBed()) return; + 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(); - if(Ke) Ke->AddToGlobal(Kgg,NULL,NULL); - delete Ke; delete tria->matice; delete tria; + + return Ke; } /*}}}*/ /*FUNCTION Penta::CreateKMatrixCouplingPattynStokes{{{1*/ -void Penta::CreateKMatrixCouplingPattynStokes( Mat Kgg){ +ElementMatrix* Penta::CreateKMatrixCouplingPattynStokes(void){ int i,j; @@ -2233,5 +2200,7 @@ if(Ke) Ke->AddToGlobal(Kgg,NULL,NULL); delete Ke; - CreateKMatrixCouplingMacAyealPattyn(Kgg); + Ke=CreateKMatrixCouplingMacAyealPattyn(); + if(Ke) Ke->AddToGlobal(Kgg,NULL,NULL); + delete Ke; break; case PattynStokesApproximationEnum: @@ -2242,5 +2211,7 @@ if(Ke) Ke->AddToGlobal(Kgg,NULL,NULL); delete Ke; - CreateKMatrixCouplingPattynStokes( Kgg); + Ke=CreateKMatrixCouplingPattynStokes(); + if(Ke) Ke->AddToGlobal(Kgg,NULL,NULL); + delete Ke; break; default: Index: /issm/trunk/src/c/objects/Elements/Penta.h =================================================================== --- /issm/trunk/src/c/objects/Elements/Penta.h (revision 5871) +++ /issm/trunk/src/c/objects/Elements/Penta.h (revision 5872) @@ -123,8 +123,8 @@ ElementMatrix* CreateKMatrixBalancedthickness(void); ElementMatrix* CreateKMatrixBalancedvelocities(void); - void CreateKMatrixCouplingMacAyealPattyn( Mat Kgg); - void CreateKMatrixCouplingMacAyealPattynViscous( Mat Kgg); - void CreateKMatrixCouplingMacAyealPattynFriction( Mat Kgg); - void CreateKMatrixCouplingPattynStokes( Mat Kgg); + ElementMatrix* CreateKMatrixCouplingMacAyealPattyn(void); + ElementMatrix* CreateKMatrixCouplingMacAyealPattynViscous(void); + ElementMatrix* CreateKMatrixCouplingMacAyealPattynFriction(void); + ElementMatrix* CreateKMatrixCouplingPattynStokes(void); void CreateKMatrixDiagnosticHoriz( Mat Kgg); ElementMatrix* CreateKMatrixDiagnosticHutter(void); Index: /issm/trunk/src/c/objects/Elements/Tria.cpp =================================================================== --- /issm/trunk/src/c/objects/Elements/Tria.cpp (revision 5871) +++ /issm/trunk/src/c/objects/Elements/Tria.cpp (revision 5872) @@ -2943,5 +2943,4 @@ } - for(i=0;ivalues[i*numdoftotal+(numdof+j)]+=Ke_gg[i][j]; for(i=0;ivalues[(i+numdof)*numdoftotal+j]+=Ke_gg[i][j];