Index: /issm/trunk/src/c/objects/Elements/Penta.cpp =================================================================== --- /issm/trunk/src/c/objects/Elements/Penta.cpp (revision 9216) +++ /issm/trunk/src/c/objects/Elements/Penta.cpp (revision 9217) @@ -783,5 +783,5 @@ ElementMatrix* Ke1=new ElementMatrix(pentabase->nodes,NUMVERTICES,this->parameters,MacAyealApproximationEnum); ElementMatrix* Ke2=new ElementMatrix(this->nodes ,NUMVERTICES,this->parameters,PattynApproximationEnum); - ElementMatrix* Ke=new ElementMatrix(Ke1,Ke2); + ElementMatrix* Ke =new ElementMatrix(Ke1,Ke2); delete Ke1; delete Ke2; @@ -1451,5 +1451,4 @@ double B[5][numdof]; double Bprime[5][numdof]; - double Ke_gg_gaussian[numdof][numdof]; //stiffness matrix evaluated at the gaussian point. Tria* tria=NULL; GaussPenta *gauss=NULL; @@ -1489,7 +1488,5 @@ &D[0][0],5,5,0, &Bprime[0][0],5,numdof,0, - &Ke_gg_gaussian[0][0],0); - - for(i=0;ivalues[i*numdof+j]+=Ke_gg_gaussian[i][j]; + &Ke->values[0],1); } @@ -1517,5 +1514,4 @@ double DL[2][2]={{ 0,0 },{0,0}}; //for basal drag double DL_scalar; - double Ke_gg_gaussian[numdof][numdof]; //stiffness matrix contribution from drag Friction *friction = NULL; GaussPenta *gauss=NULL; @@ -1565,7 +1561,5 @@ &DL[0][0],2,2,0, &L[0][0],2,numdof,0, - &Ke_gg_gaussian[0][0],0); - - for(i=0;ivalues[i*numdof+j]+=Ke_gg_gaussian[i][j]; + &Ke->values[0],1); } @@ -1780,5 +1774,4 @@ double Bprime[NDOF1][NUMVERTICES]; double DL_scalar; - double Ke_gg[numdof][numdof]={0.0}; GaussPenta *gauss=NULL; @@ -1804,7 +1797,5 @@ &DL_scalar,1,1,0, &Bprime[0][0],1,NUMVERTICES,0, - &Ke_gg[0][0],0); - - for(i=0;ivalues[i*numdof+j]+=Ke_gg[i][j]; + &Ke->values[0],1); } @@ -1830,5 +1821,4 @@ double basis[NUMVERTICES]; GaussPenta *gauss=NULL; - double Ke_g[numdof][numdof]; //stiffness matrix evaluated at the gaussian point. /*Initialize Element matrix*/ @@ -1849,14 +1839,10 @@ GetNodalFunctionsP1(&basis[0], gauss); - /**********************Do not forget the sign**********************************/ - DL_scalar=- gauss->weight*Jdet2d*surface_normal[2]; - /******************************************************************************/ + DL_scalar= - gauss->weight*Jdet2d*surface_normal[2]; TripleMultiply( basis,1,numdof,1, &DL_scalar,1,1,0, basis,1,numdof,0, - &Ke_g[0][0],0); - - for(i=0;ivalues[i*numdof+j]+=Ke_g[i][j]; + &Ke->values[0],1); } @@ -1908,8 +1894,4 @@ double D[3][3]; double K[2][2]={0.0}; - double Ke_gaussian_conduct[numdof][numdof]; - double Ke_gaussian_advec[numdof][numdof]; - double Ke_gaussian_artdiff[numdof][numdof]; - double Ke_gaussian_transient[numdof][numdof]; Tria* tria=NULL; GaussPenta *gauss=NULL; @@ -1965,5 +1947,5 @@ &D[0][0],3,3,0, &B_conduct[0][0],3,numdof,0, - &Ke_gaussian_conduct[0][0],0); + &Ke->values[0],1); /*Advection: */ @@ -1989,5 +1971,5 @@ &D[0][0],3,3,0, &Bprime_advec[0][0],3,numdof,0, - &Ke_gaussian_advec[0][0],0); + &Ke->values[0],1); /*Transient: */ @@ -2001,8 +1983,5 @@ &D_scalar_trans,1,1,0, &L[0],1,numdof,0, - &Ke_gaussian_transient[0][0],0); - } - else{ - for(i=0;ivalues[0],1); } @@ -2021,5 +2000,5 @@ &K[0][0],2,2,0, &B_artdiff[0][0],2,numdof,0, - &Ke_gaussian_artdiff[0][0],0); + &Ke->values[0],1); } else if(artdiff==2){ @@ -2031,5 +2010,5 @@ for(i=0;ivalues[i*numdof+j]=tau_parameter*D_scalar_advec* ((u-um)*dbasis[0][i]+(v-vm)*dbasis[1][i]+(w-wm)*dbasis[2][i])*((u-um)*dbasis[0][j]+(v-vm)*dbasis[1][j]+(w-wm)*dbasis[2][j]); } @@ -2038,14 +2017,9 @@ for(i=0;ivalues[i*numdof+j]+=tau_parameter*D_scalar_trans*L[j]*((u-um)*dbasis[0][i]+(v-vm)*dbasis[1][i]+(w-wm)*dbasis[2][i]); } } } } - else{ - for(i=0;ivalues[i*numdof+j]+=Ke_gaussian_conduct[i][j]+Ke_gaussian_advec[i][j]+Ke_gaussian_transient[i][j]+Ke_gaussian_artdiff[i][j]; } @@ -2070,5 +2044,4 @@ double basis[NUMVERTICES]; double D_scalar; - double Ke_gaussian[numdof][numdof]={0.0}; GaussPenta *gauss=NULL; @@ -2102,7 +2075,5 @@ &D_scalar,1,1,0, &basis[0],1,numdof,0, - &Ke_gaussian[0][0],0); - - for(i=0;ivalues[i*numdof+j]+=Ke_gaussian[i][j]; + &Ke->values[0],1); } @@ -2198,8 +2169,4 @@ double D[3][3]; double K[2][2]={0.0}; - double Ke_gaussian_conduct[numdof][numdof]; - double Ke_gaussian_advec[numdof][numdof]; - double Ke_gaussian_artdiff[numdof][numdof]; - double Ke_gaussian_transient[numdof][numdof]; Tria* tria=NULL; GaussPenta *gauss=NULL; @@ -2248,5 +2215,5 @@ &D[0][0],3,3,0, &B_conduct[0][0],3,numdof,0, - &Ke_gaussian_conduct[0][0],0); + &Ke->values[0],1); /*Advection: */ @@ -2272,5 +2239,5 @@ &D[0][0],3,3,0, &Bprime_advec[0][0],3,numdof,0, - &Ke_gaussian_advec[0][0],0); + &Ke->values[0],1); /*Transient: */ @@ -2284,8 +2251,5 @@ &D_scalar_trans,1,1,0, &L[0],1,numdof,0, - &Ke_gaussian_transient[0][0],0); - } - else{ - for(i=0;ivalues[0],1); } @@ -2304,5 +2268,5 @@ &K[0][0],2,2,0, &B_artdiff[0][0],2,numdof,0, - &Ke_gaussian_artdiff[0][0],0); + &Ke->values[0],1); } else if(artdiff==2){ @@ -2314,5 +2278,5 @@ for(i=0;ivalues[i*numdof+j]=tau_parameter*D_scalar_advec* ((u-um)*dbasis[0][i]+(v-vm)*dbasis[1][i]+(w-wm)*dbasis[2][i])*((u-um)*dbasis[0][j]+(v-vm)*dbasis[1][j]+(w-wm)*dbasis[2][j]); } @@ -2321,14 +2285,9 @@ for(i=0;ivalues[i*numdof+j]+=tau_parameter*D_scalar_trans*L[j]*((u-um)*dbasis[0][i]+(v-vm)*dbasis[1][i]+(w-wm)*dbasis[2][i]); } } } } - else{ - for(i=0;ivalues[i*numdof+j]+=Ke_gaussian_conduct[i][j]+Ke_gaussian_advec[i][j]+Ke_gaussian_transient[i][j]+Ke_gaussian_artdiff[i][j]; } @@ -2354,5 +2313,4 @@ double basis[NUMVERTICES]; double D_scalar; - double Ke_gaussian[numdof][numdof]={0.0}; GaussPenta *gauss=NULL; @@ -2386,7 +2344,5 @@ &D_scalar,1,1,0, &basis[0],1,numdof,0, - &Ke_gaussian[0][0],0); - - for(i=0;ivalues[i*numdof+j]+=Ke_gaussian[i][j]; + &Ke->values[0],1); } Index: /issm/trunk/src/c/objects/Elements/Tria.cpp =================================================================== --- /issm/trunk/src/c/objects/Elements/Tria.cpp (revision 9216) +++ /issm/trunk/src/c/objects/Elements/Tria.cpp (revision 9217) @@ -484,7 +484,4 @@ double DLprime[2][2] = {0.0}; double DL_scalar; - double Ke_gg_gaussian[numdof][numdof] = {0.0}; - double Ke_gg_thickness1[numdof][numdof] = {0.0}; - double Ke_gg_thickness2[numdof][numdof] = {0.0}; GaussTria *gauss = NULL; @@ -548,13 +545,10 @@ &DL[0][0],2,2,0, &B[0][0],2,numdof,0, - &Ke_gg_thickness1[0][0],0); + &Ke->values[0],1); TripleMultiply( &B[0][0],2,numdof,1, &DLprime[0][0],2,2,0, &Bprime[0][0],2,numdof,0, - &Ke_gg_thickness2[0][0],0); - - for(i=0;ivalues[i*numdof+j]+=Ke_gg_thickness1[i][j]; - for(i=0;ivalues[i*numdof+j]+=Ke_gg_thickness2[i][j]; + &Ke->values[0],1); if(artdiff){ @@ -565,7 +559,5 @@ &KDL[0][0],2,2,0, &Bprime[0][0],2,numdof,0, - &Ke_gg_gaussian[0][0],0); - - for(i=0;ivalues[i*numdof+j]+=Ke_gg_gaussian[i][j]; + &Ke->values[0],1); } } @@ -590,5 +582,4 @@ double DL[2][2]={0.0}; double DL_scalar; - double Ke_gg[numdof][numdof]={0.0}; GaussTria *gauss=NULL; @@ -623,7 +614,5 @@ &DL[0][0],2,2,0, &Bprime[0][0],2,numdof,0, - &Ke_gg[0][0],0); - - for(i=0;ivalues[i*numdof+j]+=Ke_gg[i][j]; + &Ke->values[0],1); } @@ -655,6 +644,4 @@ double DLprime[2][2]={0.0}; double DL_scalar; - double Ke_gg_gaussian[numdof][numdof] = {0.0}; - double Ke_gg_velocities[numdof][numdof] = {0.0}; GaussTria *gauss=NULL; @@ -734,7 +721,5 @@ &DLprime[0][0],2,2,0, &Bprime[0][0],2,numdof,0, - &Ke_gg_velocities[0][0],0); - - for(i=0;ivalues[i*numdof+j]+=Ke_gg_velocities[i][j]; + &Ke->values[0],1); if(artdiff){ @@ -745,7 +730,5 @@ &KDL[0][0],2,2,0, &Bprime[0][0],2,numdof,0, - &Ke_gg_gaussian[0][0],0); - - for(i=0;ivalues[i*numdof+j]+=Ke_gg_gaussian[i][j]; + &Ke->values[0],1); } } @@ -786,5 +769,4 @@ double D[3][3] = {0.0}; double D_scalar; - double Ke_g[numdof][numdof]; GaussTria *gauss = NULL; @@ -824,7 +806,5 @@ &D[0][0],3,3,0, &Bprime[0][0],3,numdof,0, - &Ke_g[0][0],0); - - for(i=0;ivalues[i*numdof+j]+=Ke_g[i][j]; + &Ke->values[0],1); } @@ -849,5 +829,4 @@ double L[2][numdof]; double DL[2][2] = {{ 0,0 },{0,0}}; - double Ke_g[numdof][numdof]; double DL_scalar; double slope[2] = {0.0,0.0}; @@ -894,7 +873,5 @@ &DL[0][0],2,2,0, &L[0][0],2,numdof,0, - &Ke_g[0][0],0); - - for(i=0;ivalues[i*numdof+j]+=Ke_g[i][j]; + &Ke->values[0],1); } @@ -985,7 +962,4 @@ double dLy[NUMVERTICES]; double K[2]; - double Ke_gg_gaussian1[numdof][numdof] ={0.0}; - double Ke_gg_gaussian2[numdof][numdof] ={0.0}; - double Ke_gg_gaussian3[numdof][numdof] ={0.0}; double hydro_p,hydro_q,hydro_CR,hydro_n,rho_water,g,gamma; double DL_scalar1; @@ -1019,17 +993,9 @@ GetNodalFunctions(&L[0],gauss); GetNodalFunctionsDerivatives(&dL[0][0],&xyz_list[0][0],gauss); - //printf("NUMVERTICES=%d, Jdettria=%f\n",NUMVERTICES,Jdettria); //BK for(i=0;ivalues[0],1); } if(dt){ - DL_scalar2=-gauss->weight*Jdettria*gamma*K[0]*dt; //don't forget the - - DL_scalar3=-gauss->weight*Jdettria*gamma*K[1]*dt; //don't forget the - + DL_scalar2=-gauss->weight*Jdettria*gamma*K[0]*dt; + DL_scalar3=-gauss->weight*Jdettria*gamma*K[1]*dt; } else{ - DL_scalar2=-gauss->weight*Jdettria*gamma*K[0]; //don't forget the - - DL_scalar3=-gauss->weight*Jdettria*gamma*K[1]; //don't forget the - + DL_scalar2=-gauss->weight*Jdettria*gamma*K[0]; + DL_scalar3=-gauss->weight*Jdettria*gamma*K[1]; } - TripleMultiply( &dLx[0],1,numdof,1, &DL_scalar2,1,1,0, &L[0],1,numdof,0, &Ke_gg_gaussian2[0][0],0); - TripleMultiply( &dLy[0],1,numdof,1, &DL_scalar3,1,1,0, &L[0],1,numdof,0, &Ke_gg_gaussian3[0][0],0); - - for(i=0;ivalues[i*numdof+j]+=Ke_gg_gaussian1[i][j]; - for(i=0;ivalues[i*numdof+j]+=Ke_gg_gaussian2[i][j]; - for(i=0;ivalues[i*numdof+j]+=Ke_gg_gaussian3[i][j]; - + TripleMultiply(&dLx[0],1,numdof,1, &DL_scalar2,1,1,0, &L[0],1,numdof,0,&Ke->values[0],1); + TripleMultiply(&dLy[0],1,numdof,1, &DL_scalar3,1,1,0, &L[0],1,numdof,0,&Ke->values[0],1); } @@ -1210,6 +1171,4 @@ double DLprime[2][2]={0.0}; double DL_scalar; - double Ke_gg1[numdof][numdof]={0.0}; - double Ke_gg2[numdof][numdof]={0.0}; GaussTria *gauss=NULL; @@ -1249,5 +1208,5 @@ &DL_scalar,1,1,0, &L[0],1,numdof,0, - &Ke_gg1[0][0],0); + &Ke->values[0],1); /*WARNING: B and Bprime are inverted compared to usual prognostic!!!!*/ @@ -1263,8 +1222,5 @@ &DLprime[0][0],2,2,0, &Bprime[0][0],2,numdof,0, - &Ke_gg2[0][0],0); - - for(i=0;ivalues[i*numdof+j]+=Ke_gg1[i][j]; - for(i=0;ivalues[i*numdof+j]+=Ke_gg2[i][j]; + &Ke->values[0],1); } @@ -1285,5 +1241,4 @@ double xyz_list[NUMVERTICES][3]; double L[1][3]; - double Ke_g[numdof][numdof]; GaussTria *gauss = NULL; @@ -1307,7 +1262,5 @@ &DL_scalar,1,1,0, &L[0][0],1,3,0, - &Ke_g[0][0],0); - - for(i=0;ivalues[i*numdof+j]+=Ke_g[i][j]; + &Ke->values[0],1); } @@ -2868,7 +2821,4 @@ K[1]=pow(w,hydro_p)*pow((rho_ice*g*dsdy+(rho_water/rho_ice-1)*rho_ice*g*dbdy) - rho_ice * g * kn/w * (dsdy - dbdy ) * surface_slope,hydro_q); - //printf("K[0]=%g,K[1]=%g,w=%g\n",K[0],K[1],w); //bk - //if (w<0) {printf("negative w!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!\n"); - //printf("K[0]=%g,K[1]=%g,w=%g\n",K[0],K[1],w);} } /*}}}*/