Index: issm/trunk/src/c/objects/Elements/Penta.cpp =================================================================== --- issm/trunk/src/c/objects/Elements/Penta.cpp (revision 5848) +++ issm/trunk/src/c/objects/Elements/Penta.cpp (revision 5849) @@ -2200,5 +2200,7 @@ break; case PattynApproximationEnum: - CreateKMatrixDiagnosticPattyn( Kgg); + Ke=CreateKMatrixDiagnosticPattyn(); + if(Ke) Ke->AddToGlobal(Kgg,NULL,NULL); + delete Ke; break; case StokesApproximationEnum: @@ -2211,9 +2213,13 @@ case MacAyealPattynApproximationEnum: CreateKMatrixDiagnosticMacAyeal3d( Kgg); - CreateKMatrixDiagnosticPattyn( Kgg); - CreateKMatrixCouplingMacAyealPattyn( Kgg); + Ke=CreateKMatrixDiagnosticPattyn(); + if(Ke) Ke->AddToGlobal(Kgg,NULL,NULL); + delete Ke; + CreateKMatrixCouplingMacAyealPattyn(Kgg); break; case PattynStokesApproximationEnum: - CreateKMatrixDiagnosticPattyn( Kgg); + Ke=CreateKMatrixDiagnosticPattyn(); + if(Ke) Ke->AddToGlobal(Kgg,NULL,NULL); + delete Ke; CreateKMatrixDiagnosticStokes( Kgg); CreateKMatrixCouplingPattynStokes( Kgg); @@ -2415,28 +2421,31 @@ /*}}}*/ /*FUNCTION Penta::CreateKMatrixDiagnosticPattyn{{{1*/ -void Penta::CreateKMatrixDiagnosticPattyn( Mat Kgg){ +ElementMatrix* Penta::CreateKMatrixDiagnosticPattyn(void){ + + /*compute all stiffness matrices for this element*/ + ElementMatrix* Ke1=CreateKMatrixDiagnosticPattynViscous(); + ElementMatrix* Ke2=CreateKMatrixDiagnosticPattynFriction(); + ElementMatrix* Ke =new ElementMatrix(Ke1,Ke2); + + /*clean-up and return*/ + delete Ke1; + delete Ke2; + return Ke; +} +/*}}}*/ +/*FUNCTION Penta::CreateKMatrixDiagnosticPattynViscous{{{1*/ +ElementMatrix* Penta::CreateKMatrixDiagnosticPattynViscous(void){ /* local declarations */ - int i,j; - - /* node data: */ const int numdof=2*NUMVERTICES; + int i,j,ig; double xyz_list[NUMVERTICES][3]; - int* doflist=NULL; - - /* 3d gaussian points: */ - int ig; GaussPenta *gauss=NULL; - - /* material data: */ double viscosity; //viscosity double oldviscosity; //viscosity double newviscosity; //viscosity - - /* strain rate: */ + double viscosity_overshoot; double epsilon[5]; /* epsilon=[exx,eyy,exy,exz,eyz];*/ double oldepsilon[5]; /* epsilon=[exx,eyy,exy,exz,eyz];*/ - - /* matrices: */ double B[5][numdof]; double Bprime[5][numdof]; @@ -2446,40 +2455,22 @@ double DL[2][2]={0.0}; //for basal drag double DL_scalar; - - /* local element matrices: */ 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 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; - /*retrieve some parameters: */ + /*Initialize Element matrix and return if necessary*/ + if(IsOnWater()) return NULL; + ElementMatrix* Ke=this->NewElementMatrix(PattynApproximationEnum); + + /*Retrieve all inputs and parameters*/ + GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES); this->parameters->FindParam(&viscosity_overshoot,ViscosityOvershootEnum); - - /*If on water, skip stiffness: */ - if(IsOnWater())return; - - /*Implement standard penta element: */ - - /* Get node coordinates and dof list: */ - GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES); - GetDofList(&doflist,PattynApproximationEnum,GsetEnum); - - /*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); @@ -2493,21 +2484,13 @@ gauss->GaussPoint(ig); - /*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: */ GetBPattyn(&B[0][0], &xyz_list[0][0], gauss); GetBprimePattyn(&Bprime[0][0], &xyz_list[0][0], gauss); - /* 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); @@ -2515,5 +2498,4 @@ for (i=0;i<5;i++) D[i][i]=D_scalar; - /* Do the triple product tB*D*Bprime: */ TripleMultiply( &B[0][0],5,numdof,1, &D[0][0],5,5,0, @@ -2521,26 +2503,27 @@ &Ke_gg_gaussian[0][0],0); - /* Add the Ke_gg_gaussian, and optionally Ke_gg_gaussian onto Ke_gg: */ - for(i=0;iCreateKMatrixDiagnosticPattynFriction(Kgg); - delete tria->matice; delete tria; + for(i=0;ivalues[i*numdof+j]+=Ke_gg_gaussian[i][j]; } /*Clean up and return*/ delete gauss; - xfree((void**)&doflist); + return Ke; +} +/*}}}*/ +/*FUNCTION Penta::CreateKMatrixDiagnosticPattynFriction{{{1*/ +ElementMatrix* Penta::CreateKMatrixDiagnosticPattynFriction(void){ + + /*Initialize Element matrix and return if necessary*/ + if(IsOnWater() || IsOnShelf() || !IsOnBed()) return NULL; + + /*Build a tria element using the 3 grids of the base of the penta. Then use + * the tria functionality to build a friction stiffness matrix on these 3 + * grids: */ + Tria* tria=(Tria*)SpawnTria(0,1,2); //grids 0, 1 and 2 make the new tria. + ElementMatrix* Ke=tria->CreateKMatrixDiagnosticPattynFriction(); + delete tria->matice; delete tria; + + /*clean-up and return*/ + return Ke; } /*}}}*/ @@ -4084,5 +4067,5 @@ /*First, figure out size of doflist: */ numdof=0; - for(i=0;i<3;i++){ + for(i=0;iGetNumberOfDofs(approximation_enum,setenum); numdof+=ndof_list[i]; @@ -4095,5 +4078,5 @@ /*Populate: */ count=0; - for(i=0;i<3;i++){ + for(i=0;iGetDofList(&doflist[count],approximation_enum,setenum); count+=ndof_list[i]; Index: issm/trunk/src/c/objects/Elements/Penta.h =================================================================== --- issm/trunk/src/c/objects/Elements/Penta.h (revision 5848) +++ issm/trunk/src/c/objects/Elements/Penta.h (revision 5849) @@ -129,5 +129,7 @@ ElementMatrix* CreateKMatrixDiagnosticMacAyeal2d(void); void CreateKMatrixDiagnosticMacAyeal3d(Mat Kgg); - void CreateKMatrixDiagnosticPattyn( Mat Kgg); + ElementMatrix* CreateKMatrixDiagnosticPattyn(void); + ElementMatrix* CreateKMatrixDiagnosticPattynViscous(void); + ElementMatrix* CreateKMatrixDiagnosticPattynFriction(void); void CreateKMatrixDiagnosticStokes( Mat Kgg); void CreateKMatrixDiagnosticVert( Mat Kgg); Index: issm/trunk/src/c/objects/Elements/Tria.cpp =================================================================== --- issm/trunk/src/c/objects/Elements/Tria.cpp (revision 5848) +++ issm/trunk/src/c/objects/Elements/Tria.cpp (revision 5849) @@ -2978,49 +2978,34 @@ /*}}}*/ /*FUNCTION Tria::CreateKMatrixDiagnosticPattynFriction {{{1*/ -void Tria::CreateKMatrixDiagnosticPattynFriction(Mat Kgg){ - - /* local declarations */ +ElementMatrix* Tria::CreateKMatrixDiagnosticPattynFriction(void){ + + const int numdof = NDOF2*NUMVERTICES; int i,j; + int ig; int analysis_type; - - /* node data: */ - const int numdof = 2 *NUMVERTICES; double xyz_list[NUMVERTICES][3]; - int* doflist=NULL; int numberofdofspernode=2; - - /* gaussian points: */ - int ig; GaussTria *gauss=NULL; - - /* matrices: */ double L[2][numdof]; double DL[2][2]={{ 0,0 },{0,0}}; //for basal drag double DL_scalar; - - /* local element matrices: */ double Ke_gg[numdof][numdof]={0.0}; double Ke_gg_gaussian[numdof][numdof]; //stiffness matrix contribution from drag - double Jdet; - - /*slope: */ double slope[2]={0.0,0.0}; double slope_magnitude; - - /*friction: */ Friction *friction = NULL; double alpha2; - double MAXSLOPE=.06; // 6 % double MOUNTAINKEXPONENT=10; - - /*inputs: */ int drag_type; - /*retrive parameters: */ + /*Initialize Element matrix and return if necessary*/ + if(IsOnWater() || IsOnShelf()) return NULL; + ElementMatrix* Ke=this->NewElementMatrix(PattynApproximationEnum); + + /*Retrieve all inputs and parameters*/ + GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES); parameters->FindParam(&analysis_type,AnalysisTypeEnum); - - /*retrieve inputs :*/ inputs->GetParameterValue(&drag_type,DragTypeEnum); Input* surface_input=inputs->GetInput(SurfaceEnum); ISSMASSERT(surface_input); @@ -3028,13 +3013,4 @@ Input* vy_input=inputs->GetInput(VyEnum); ISSMASSERT(vy_input); Input* vz_input=inputs->GetInput(VzEnum); ISSMASSERT(vz_input); - - /* Get node coordinates and dof list: */ - GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES); - GetDofList(&doflist,PattynApproximationEnum,GsetEnum); - - if (IsOnShelf()){ - /*no friction, do nothing*/ - return; - } /*build friction object, used later on: */ @@ -3048,5 +3024,4 @@ gauss->GaussPoint(ig); - /*Friction: */ friction->GetAlpha2(&alpha2, gauss,VxEnum,VyEnum,VzEnum); // TO UNCOMMENT @@ -3060,17 +3035,10 @@ } - /* Get Jacobian determinant: */ + GetL(&L[0][0], &xyz_list[0][0], gauss,numberofdofspernode); GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss); - - /*Get L matrix: */ - GetL(&L[0][0], &xyz_list[0][0], gauss,numberofdofspernode); - DL_scalar=alpha2*gauss->weight*Jdet; - for (i=0;i<2;i++){ - DL[i][i]=DL_scalar; - } + for (i=0;i<2;i++) DL[i][i]=DL_scalar; - /* Do the triple producte tL*D*L: */ TripleMultiply( &L[0][0],2,numdof,1, &DL[0][0],2,2,0, @@ -3078,15 +3046,11 @@ &Ke_gg_gaussian[0][0],0); - for( i=0; ivalues[i*numdof+j]+=Ke_gg_gaussian[i][j]; + } /*Clean up and return*/ delete gauss; delete friction; - xfree((void**)&doflist); + return Ke; } /*}}}*/ @@ -5118,5 +5082,5 @@ /*First, figure out size of doflist: */ numdof=0; - for(i=0;i<3;i++){ + for(i=0;iGetNumberOfDofs(approximation_enum,setenum); numdof+=ndof_list[i]; @@ -5129,5 +5093,5 @@ /*Populate: */ count=0; - for(i=0;i<3;i++){ + for(i=0;iGetDofList(&doflist[count],approximation_enum,setenum); count+=ndof_list[i]; Index: issm/trunk/src/c/objects/Elements/Tria.h =================================================================== --- issm/trunk/src/c/objects/Elements/Tria.h (revision 5848) +++ issm/trunk/src/c/objects/Elements/Tria.h (revision 5849) @@ -130,5 +130,5 @@ void CreateKMatrixDiagnosticMacAyealFriction(Mat Kgg,Mat Kff, Mat Kfs); void CreateKMatrixCouplingMacAyealPattynFriction(Mat Kgg); - void CreateKMatrixDiagnosticPattynFriction(Mat Kgg); + ElementMatrix* CreateKMatrixDiagnosticPattynFriction(void); ElementMatrix* CreateKMatrixDiagnosticHutter(void); void CreateKMatrixDiagnosticSurfaceVert(Mat Kgg); Index: issm/trunk/src/c/shared/Elements/Paterson.cpp =================================================================== --- issm/trunk/src/c/shared/Elements/Paterson.cpp (revision 5848) +++ issm/trunk/src/c/shared/Elements/Paterson.cpp (revision 5849) @@ -32,32 +32,32 @@ B=pow((double)10,(double)8)*(-0.000292866376675*pow(T+50,3)+ 0.011672640664130*pow(T+50,2) -0.325004442485481*(T+50)+ 6.524779401948101); } - if((T>=-45.0) && (T<=-40.0)){ + else if((T>=-45.0) && (T<=-40.0)){ B=pow((double)10,(double)8)*(-0.000292866376675*pow(T+45,3)+ 0.007279645014004*pow(T+45,2) -0.230243014094813*(T+45)+ 5.154964909039554); } - if((T>=-40.0) && (T<=-35.0)){ + else if((T>=-40.0) && (T<=-35.0)){ B=pow((double)10,(double)8)*(0.000072737147457*pow(T+40,3)+ 0.002886649363879*pow(T+40,2) -0.179411542205399*(T+40)+ 4.149132666831214); } - if((T>=-35.0) && (T<=-30.0)){ + else if((T>=-35.0) && (T<=-30.0)){ B=pow((double)10,(double)8)*(-0.000086144770023*pow(T+35,3)+ 0.003977706575736*pow(T+35,2) -0.145089762507325*(T+35)+ 3.333333333333331); } - if((T>=-30.0) && (T<=-25.0)){ + else if((T>=-30.0) && (T<=-25.0)){ B=pow((double)10,(double)8)*(-0.000043984685769*pow(T+30,3)+ 0.002685535025386*pow(T+30,2) -0.111773554501713*(T+30)+ 2.696559088937191); } - if((T>=-25.0) && (T<=-20.0)){ + else if((T>=-25.0) && (T<=-20.0)){ B=pow((double)10,(double)8)*(-0.000029799523463*pow(T+25,3)+ 0.002025764738854*pow(T+25,2) -0.088217055680511*(T+25)+ 2.199331606342181); } - if((T>=-20.0) && (T<=-15.0)){ + else if((T>=-20.0) && (T<=-15.0)){ B=pow((double)10,(double)8)*(0.000136920904777*pow(T+20,3)+ 0.001578771886910*pow(T+20,2) -0.070194372551690*(T+20)+ 1.805165505978111); } - if((T>=-15.0) && (T<=-10.0)){ + else if((T>=-15.0) && (T<=-10.0)){ B=pow((double)10,(double)8)*(-0.000899763781026*pow(T+15,3)+ 0.003632585458564*pow(T+15,2) -0.044137585824322*(T+15)+ 1.510778053489523); } - if((T>=-10.0) && (T<=-5.0)){ + else if((T>=-10.0) && (T<=-5.0)){ B=pow((double)10,(double)8)*(0.001676964325070*pow(T+10,3)- 0.009863871256831*pow(T+10,2) -0.075294014815659*(T+10)+ 1.268434288203714); } - if((T>=-5.0) && (T<=-2.0)){ + else if((T>=-5.0) && (T<=-2.0)){ B=pow((double)10,(double)8)*(-0.003748937622487*pow(T+5,3)+0.015290593619213*pow(T+5,2) -0.048160403003748*(T+5)+ 0.854987973338348); } - if(T>=-2.0){ + else if(T>=-2.0){ B=pow((double)10,(double)8)*(-0.003748937622488*pow(T+2,3)-0.018449844983174*pow(T+2,2) -0.057638157095631*(T+2)+ 0.746900791092860); }