Index: /issm/trunk/src/c/objects/Elements/Tria.cpp =================================================================== --- /issm/trunk/src/c/objects/Elements/Tria.cpp (revision 5843) +++ /issm/trunk/src/c/objects/Elements/Tria.cpp (revision 5844) @@ -2422,5 +2422,4 @@ double DLprime[2][2] = {0.0}; double DL_scalar; - double Ke_gg[numdof][numdof] = {0.0}; double Ke_gg_gaussian[numdof][numdof] = {0.0}; double Ke_gg_thickness1[numdof][numdof] = {0.0}; @@ -2428,5 +2427,5 @@ GaussTria *gauss = NULL; - /*Initialize Element matrix and return if necessary*/ + /*Initialize Element matrix and return if necessary*/ if(IsOnWater()) return NULL; ElementMatrix* Ke=this->NewElementMatrix(NoneApproximationEnum); @@ -2447,5 +2446,5 @@ } - //Create Artificial diffusivity once for all if requested + /*Create Artificial diffusivity once for all if requested*/ if(artdiff){ gauss=new GaussTria(); @@ -2460,5 +2459,5 @@ } - /* Start looping on the number of gaussian points: */ + /*Start looping on the number of gaussian points:*/ gauss=new GaussTria(2); for (ig=gauss->begin();igend();ig++){ @@ -2519,20 +2518,17 @@ ElementMatrix* Tria::CreateKMatrixBalancedthickness_DG(void){ - /* local declarations */ - const int numdof=NDOF1*NUMVERTICES; - int i,j; - int ig; - double xyz_list[NUMVERTICES][3]; - GaussTria *gauss=NULL; - double B[2][NUMVERTICES]; - double Bprime[2][NUMVERTICES]; - double DL[2][2]={0.0}; - double DLprime[2][2]={0.0}; - double DL_scalar; - double Ke_gg[numdof][numdof]={0.0}; - double Ke_gg2[numdof][numdof]={0.0}; - double Jdettria; - double vx,vy; - int dim; + /*Constants*/ + const int numdof=NDOF1*NUMVERTICES; + + /*Intermediaries*/ + int i,j,ig,dim; + double vx,vy,Jdettria; + double xyz_list[NUMVERTICES][3]; + double B[2][NUMVERTICES]; + double Bprime[2][NUMVERTICES]; + double DL[2][2]={0.0}; + double DL_scalar; + double Ke_gg[numdof][numdof]={0.0}; + GaussTria *gauss=NULL; /*Initialize Element matrix and return if necessary*/ @@ -2546,5 +2542,5 @@ Input* vy_input=inputs->GetInput(VyEnum); ISSMASSERT(vy_input); - /* Start looping on the number of gaussian points: */ + /*Start looping on the number of gaussian points:*/ gauss=new GaussTria(2); for (ig=gauss->begin();igend();ig++){ @@ -2552,4 +2548,5 @@ gauss->GaussPoint(ig); + GetJacobianDeterminant2d(&Jdettria, &xyz_list[0][0],gauss); /*WARNING: B and Bprime are inverted compared to usual prognostic!!!!*/ GetBPrognostic(&Bprime[0][0], &xyz_list[0][0], gauss); @@ -2559,16 +2556,14 @@ vy_input->GetParameterValue(&vy,gauss); - GetJacobianDeterminant2d(&Jdettria, &xyz_list[0][0],gauss); DL_scalar=-gauss->weight*Jdettria; - - DLprime[0][0]=DL_scalar*vx; - DLprime[1][1]=DL_scalar*vy; + DL[0][0]=DL_scalar*vx; + DL[1][1]=DL_scalar*vy; TripleMultiply( &B[0][0],2,numdof,1, - &DLprime[0][0],2,2,0, + &DL[0][0],2,2,0, &Bprime[0][0],2,numdof,0, - &Ke_gg2[0][0],0); - - for(i=0;ivalues[i*numdof+j]+=Ke_gg2[i][j]; + &Ke_gg[0][0],0); + + for(i=0;ivalues[i*numdof+j]+=Ke_gg[i][j]; } @@ -2581,33 +2576,26 @@ ElementMatrix* Tria::CreateKMatrixBalancedvelocities(void){ - /* local declarations */ + /*Constants*/ const int numdof=NDOF1*NUMVERTICES; - int i,j; - double xyz_list[NUMVERTICES][3]; - int ig; - GaussTria *gauss=NULL; - double L[NUMVERTICES]; - double B[2][NUMVERTICES]; - double Bprime[2][NUMVERTICES]; - double DL[2][2]={0.0}; - double DLprime[2][2]={0.0}; - double DL_scalar; - double Ke_gg[numdof][numdof]={0.0};//local element stiffness matrix - double Ke_gg_gaussian[numdof][numdof]={0.0}; //stiffness matrix evaluated at the gaussian point. - double Ke_gg_velocities1[numdof][numdof]={0.0}; //stiffness matrix evaluated at the gaussian point. - double Ke_gg_velocities2[numdof][numdof]={0.0}; //stiffness matrix evaluated at the gaussian point. - double Jdettria; + + /*Intermediaries */ + bool artdiff; + int i,j,ig,dim; + double nx,ny,norm,Jdettria; + double dvx[2],dvy[2]; + double vx,vy,dvxdx,dvydy; + double v_gauss[2]={0.0}; double surface_normal[3]; double surface_list[3]; - double nx,ny,norm; - double dvx[2]; - double dvy[2]; - double vx,vy; - double dvxdx,dvydy; - double v_gauss[2]={0.0}; + double xyz_list[NUMVERTICES][3]; + double B[2][NUMVERTICES]; + double Bprime[2][NUMVERTICES]; double K[2][2]={0.0}; double KDL[2][2]={0.0}; - int dim; - bool artdiff; + 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; /*Initialize Element matrix and return if necessary*/ @@ -2644,10 +2632,8 @@ } if(nx==0 && ny==0){ - nx=0; - ny=1; + nx=0; ny=1; } norm=pow( pow(nx,2)+pow(ny,2) , (double).5); - nx=nx/norm; - ny=ny/norm; + nx=nx/norm; ny=ny/norm; //Create Artificial diffusivity once for all if requested @@ -2660,5 +2646,4 @@ vxaverage_input->GetParameterAverage(&v_gauss[0]); vyaverage_input->GetParameterAverage(&v_gauss[1]); - K[0][0]=pow(10,2)*pow(Jdettria,(double).5)/2.0*fabs(v_gauss[0]); //pow should be zero!! K[1][1]=pow(10,2)*pow(Jdettria,(double).5)/2.0*fabs(v_gauss[1]); @@ -2671,7 +2656,7 @@ gauss->GaussPoint(ig); + GetJacobianDeterminant2d(&Jdettria, &xyz_list[0][0],gauss); GetBPrognostic(&B[0][0], &xyz_list[0][0], gauss); GetBprimePrognostic(&Bprime[0][0], &xyz_list[0][0], gauss); - GetJacobianDeterminant2d(&Jdettria, &xyz_list[0][0],gauss); vxaverage_input->GetParameterValue(&vx,gauss); @@ -2679,7 +2664,7 @@ vxaverage_input->GetParameterDerivativeValue(&dvx[0],&xyz_list[0][0],gauss); vyaverage_input->GetParameterDerivativeValue(&dvy[0],&xyz_list[0][0],gauss); + dvxdx=dvx[0]; dvydy=dvy[1]; - DL_scalar=gauss->weight*Jdettria; @@ -2690,7 +2675,7 @@ &DLprime[0][0],2,2,0, &Bprime[0][0],2,numdof,0, - &Ke_gg_velocities2[0][0],0); - - for(i=0;ivalues[i*numdof+j]+=Ke_gg_velocities2[i][j]; + &Ke_gg_velocities[0][0],0); + + for(i=0;ivalues[i*numdof+j]+=Ke_gg_velocities[i][j]; if(artdiff){ @@ -2704,7 +2689,5 @@ for(i=0;ivalues[i*numdof+j]+=Ke_gg_gaussian[i][j]; - } - } @@ -2738,9 +2721,7 @@ int i,j,ig; double xyz_list[NUMVERTICES][3]; - GaussTria *gauss = NULL; double viscosity,newviscosity,oldviscosity; - double viscosity_overshoot,thickness; - double epsilon[3]; /* epsilon=[exx,eyy,exy]; */ - double oldepsilon[3]; /* oldepsilon=[exx,eyy,exy]; */ + double viscosity_overshoot,thickness,Jdet; + double epsilon[3],oldepsilon[3]; /* epsilon=[exx,eyy,exy]; */ double B[3][numdof]; double Bprime[3][numdof]; @@ -2748,5 +2729,5 @@ double D_scalar; double Ke_g[numdof][numdof]; - double Jdet; + GaussTria *gauss = NULL; /*Initialize Element matrix and return if necessary*/ @@ -2769,17 +2750,17 @@ gauss->GaussPoint(ig); - thickness_input->GetParameterValue(&thickness, gauss); + GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss); + GetBMacAyeal(&B[0][0], &xyz_list[0][0], gauss); + GetBprimeMacAyeal(&Bprime[0][0], &xyz_list[0][0], gauss); + this->GetStrainRate2d(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input); this->GetStrainRate2d(&oldepsilon[0],&xyz_list[0][0],gauss,vxold_input,vyold_input); matice->GetViscosity2d(&viscosity, &epsilon[0]); matice->GetViscosity2d(&oldviscosity, &oldepsilon[0]); - GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss); + thickness_input->GetParameterValue(&thickness, gauss); newviscosity=viscosity+viscosity_overshoot*(viscosity-oldviscosity); D_scalar=2*newviscosity*thickness*gauss->weight*Jdet; for (i=0;i<3;i++) D[i][i]=D_scalar; - - GetBMacAyeal(&B[0][0], &xyz_list[0][0], gauss); - GetBprimeMacAyeal(&Bprime[0][0], &xyz_list[0][0], gauss); TripleMultiply(&B[0][0],3,numdof,1,