Index: /issm/trunk/src/c/objects/Elements/Tria.cpp =================================================================== --- /issm/trunk/src/c/objects/Elements/Tria.cpp (revision 5877) +++ /issm/trunk/src/c/objects/Elements/Tria.cpp (revision 5878) @@ -3098,15 +3098,16 @@ ElementMatrix* Tria::CreateKMatrixMelting(void){ + /*Constants*/ const int numdof=NUMVERTICES*NDOF1; - int i,j,ig; + + /*Intermediaries */ + int i,j,ig; + double heatcapacity,latentheat; + double Jdet,D_scalar; double xyz_list[NUMVERTICES][3]; - double heatcapacity,latentheat; - GaussTria *gauss=NULL; - double Jdet; - double D_scalar; - double K_terms[numdof][numdof]={0.0}; double L[3]; double tLD[3]; double Ke_gaussian[numdof][numdof]={0.0}; + GaussTria *gauss=NULL; /*Initialize Element matrix and return if necessary*/ @@ -3127,4 +3128,5 @@ GetL(&L[0], &xyz_list[0][0], gauss,NDOF1); GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0], gauss); + D_scalar=latentheat/heatcapacity*gauss->weight*Jdet; @@ -3157,32 +3159,26 @@ ElementMatrix* Tria::CreateKMatrixPrognostic_CG(void){ - /* local declarations */ + /*Constants*/ const int numdof=NDOF1*NUMVERTICES; - int i,j; - int ig; - double xyz_list[NUMVERTICES][3]; - int numberofdofspernode=1; + + /*Intermediaries */ + bool artdiff; + int i,j,ig,dim; + double Jdettria,DL_scalar,dt; + double vx,vy,dvxdx,dvydy; + double dvx[2],dvy[2]; + double v_gauss[2]={0.0}; + double xyz_list[NUMVERTICES][3]; + double L[NUMVERTICES]; + double B[2][NUMVERTICES]; + double Bprime[2][NUMVERTICES]; + double K[2][2] ={0.0}; + double KDL[2][2] ={0.0}; + double DL[2][2] ={0.0}; + double DLprime[2][2] ={0.0}; + 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; - 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}; - double Ke_gg_gaussian[numdof][numdof]={0.0}; - double Ke_gg_thickness1[numdof][numdof]={0.0}; - double Ke_gg_thickness2[numdof][numdof]={0.0}; - double Jdettria; - double dvx[2]; - double dvy[2]; - double vx,vy; - double dvxdx,dvydy; - double v_gauss[2]={0.0}; - double K[2][2]={0.0}; - double KDL[2][2]={0.0}; - int dim; - double dt; - bool artdiff; /*Initialize Element matrix and return if necessary*/ @@ -3208,5 +3204,4 @@ //Create Artificial diffusivity once for all if requested if(artdiff){ - //Get the Jacobian determinant gauss=new GaussTria(); gauss->GaussCenter(); @@ -3214,5 +3209,4 @@ delete gauss; - //Build K matrix (artificial diffusivity matrix) vxaverage_input->GetParameterAverage(&v_gauss[0]); vyaverage_input->GetParameterAverage(&v_gauss[1]); @@ -3228,6 +3222,11 @@ gauss->GaussPoint(ig); - GetL(&L[0], &xyz_list[0][0], gauss,numberofdofspernode); GetJacobianDeterminant2d(&Jdettria, &xyz_list[0][0],gauss); + GetL(&L[0], &xyz_list[0][0], gauss,NDOF1); + + vxaverage_input->GetParameterValue(&vx,gauss); + vyaverage_input->GetParameterValue(&vy,gauss); + vxaverage_input->GetParameterDerivativeValue(&dvx[0],&xyz_list[0][0],gauss); + vyaverage_input->GetParameterDerivativeValue(&dvy[0],&xyz_list[0][0],gauss); DL_scalar=gauss->weight*Jdettria; @@ -3241,8 +3240,4 @@ GetBprimePrognostic(&Bprime[0][0], &xyz_list[0][0], gauss); - vxaverage_input->GetParameterValue(&vx,gauss); - vyaverage_input->GetParameterValue(&vy,gauss); - 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]; @@ -3289,25 +3284,21 @@ ElementMatrix* Tria::CreateKMatrixPrognostic_DG(void){ - /* local declarations */ - int i,j; + /*Constants*/ const int numdof=NDOF1*NUMVERTICES; - double xyz_list[NUMVERTICES][3]; - int* doflist=NULL; - int numberofdofspernode=1; - 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}; - double Ke_gg1[numdof][numdof]={0.0}; - double Ke_gg2[numdof][numdof]={0.0}; - double Jdettria; - double vx,vy; - int dim; - double dt; + + /*Intermediaries */ + int i,j,ig,dim; + int* doflist=NULL; + double xyz_list[NUMVERTICES][3]; + double Jdettria,dt,vx,vy; + 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_gg1[numdof][numdof]={0.0}; + double Ke_gg2[numdof][numdof]={0.0}; + GaussTria *gauss=NULL; /*Initialize Element matrix and return if necessary*/ @@ -3336,7 +3327,9 @@ gauss->GaussPoint(ig); - - GetL(&L[0], &xyz_list[0][0], gauss,numberofdofspernode); + vxaverage_input->GetParameterValue(&vx,gauss); + vyaverage_input->GetParameterValue(&vy,gauss); + GetJacobianDeterminant2d(&Jdettria, &xyz_list[0][0],gauss); + GetL(&L[0], &xyz_list[0][0], gauss,NDOF1); DL_scalar=gauss->weight*Jdettria; @@ -3347,11 +3340,8 @@ &Ke_gg1[0][0],0); - /*Get B and B prime matrix: */ /*WARNING: B and Bprime are inverted compared to usual prognostic!!!!*/ GetBPrognostic(&Bprime[0][0], &xyz_list[0][0], gauss); GetBprimePrognostic(&B[0][0], &xyz_list[0][0], gauss); - vxaverage_input->GetParameterValue(&vx,gauss); - vyaverage_input->GetParameterValue(&vy,gauss); DL_scalar=-dt*gauss->weight*Jdettria; @@ -3632,27 +3622,16 @@ void Tria::CreatePVectorBalancedvelocities(Vec pg){ - /* local declarations */ - int i,j; - - /* node data: */ + /*Constants*/ const int numdof=NDOF1*NUMVERTICES; - double xyz_list[NUMVERTICES][3]; - int* doflist=NULL; - int numberofdofspernode=1; - - /* gaussian points: */ - int ig; + + /*Intermediaries */ + int i,j,ig; + int* doflist=NULL; + double xyz_list[NUMVERTICES][3]; + double Jdettria,accumulation_g,melting_g; + double L[NUMVERTICES]; + double pe_g[NUMVERTICES]={0.0}; GaussTria* gauss=NULL; - /* matrix */ - double pe_g[NUMVERTICES]={0.0}; - double L[NUMVERTICES]; - double Jdettria; - - /*input parameters for structural analysis (diagnostic): */ - double accumulation_g; - double melting_g; - - /* Get node coordinates and dof list: */ GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES); GetDofList(&doflist,NoneApproximationEnum,GsetEnum); @@ -3668,20 +3647,14 @@ gauss->GaussPoint(ig); - /* Get Jacobian determinant: */ GetJacobianDeterminant2d(&Jdettria, &xyz_list[0][0],gauss); - - /*Get L matrix: */ - GetL(&L[0], &xyz_list[0][0], gauss,numberofdofspernode); - - /* Get accumulation, melting at gauss point */ + GetL(&L[0], &xyz_list[0][0], gauss,NDOF1); + accumulation_input->GetParameterValue(&accumulation_g,gauss); melting_input->GetParameterValue(&melting_g,gauss); - /* Add value into pe_g: */ for( i=0; iweight*(accumulation_g-melting_g)*L[i]; } - /*Add pe_g to global matrix Kgg: */ VecSetValues(pg,numdof,doflist,(const double*)pe_g,ADD_VALUES);