Index: /issm/trunk/src/c/objects/Elements/Tria.cpp =================================================================== --- /issm/trunk/src/c/objects/Elements/Tria.cpp (revision 5879) +++ /issm/trunk/src/c/objects/Elements/Tria.cpp (revision 5880) @@ -3467,26 +3467,15 @@ void Tria::CreatePVectorBalancedthickness_CG(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 dhdt_g,melting_g,accumulation_g,Jdettria; + 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; - double dhdt_g; /* Get node coordinates and dof list: */ @@ -3505,21 +3494,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 */ accumulation_input->GetParameterValue(&accumulation_g,gauss); melting_input->GetParameterValue(&melting_g,gauss); dhdt_input->GetParameterValue(&dhdt_g,gauss); - /* Add value into pe_g: */ + GetJacobianDeterminant2d(&Jdettria, &xyz_list[0][0],gauss); + GetL(&L[0], &xyz_list[0][0], gauss,NDOF1); + for( i=0; iweight*(accumulation_g-melting_g-dhdt_g)*L[i]; - - } - - /*Add pe_g to global matrix Kgg: */ + } + VecSetValues(pg,numdof,doflist,(const double*)pe_g,ADD_VALUES); @@ -3532,26 +3514,15 @@ void Tria::CreatePVectorBalancedthickness_DG(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 melting_g,accumulation_g,dhdt_g,Jdettria; + 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; - double dhdt_g; /* Get node coordinates and dof list: */ @@ -3570,21 +3541,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 and thickness at gauss point */ accumulation_input->GetParameterValue(&accumulation_g,gauss); melting_input->GetParameterValue(&melting_g,gauss); dhdt_input->GetParameterValue(&dhdt_g,gauss); - /* Add value into pe_g: */ - for( i=0; iweight*(accumulation_g-melting_g-dhdt_g)*L[i]; - - } - - /*Add pe_g to global matrix Kgg: */ + GetJacobianDeterminant2d(&Jdettria, &xyz_list[0][0],gauss); + GetL(&L[0], &xyz_list[0][0], gauss,NDOF1); + + for(i=0;iweight*(accumulation_g-melting_g-dhdt_g)*L[i]; + } + VecSetValues(pg,numdof,doflist,(const double*)pe_g,ADD_VALUES); @@ -3642,34 +3606,19 @@ void Tria::CreatePVectorDiagnosticBaseVert(Vec pg){ - int i,j; - - /* node data: */ + /*Constants*/ const int numdof=NDOF1*NUMVERTICES; - double xyz_list[NUMVERTICES][3]; - int* doflist=NULL; - - /* gaussian points: */ - int ig; + + /*Intermediaries */ + int i,j,ig; + int approximation; + int* doflist=NULL; + double xyz_list[NUMVERTICES][3]; + double Jdet; + double vx,vy,vz,dbdx,dbdy,meltingvalue; + double slope[2]; + double L[NUMVERTICES]; + double pe_g[numdof]={0.0}; + double pe_g_gaussian[numdof]; GaussTria* gauss=NULL; - - /* Jacobian: */ - double Jdet; - - /*nodal functions: */ - double l1l2l3[3]; - - /*element vector at the gaussian points: */ - double pe_g[numdof]={0.0}; - double pe_g_gaussian[numdof]; - - /* matrices: */ - double L[NUMVERTICES]; - - /*input parameters for structural analysis (diagnostic): */ - double vx,vy,vz; - double meltingvalue; - double slope[2]; - double dbdx,dbdy; - int approximation; /* Get node coordinates and dof list: */ @@ -3679,4 +3628,5 @@ /*Retrieve all inputs we will be needing: */ inputs->GetParameterValue(&approximation,ApproximationEnum); + Input* bed_input=inputs->GetInput(BedEnum); ISSMASSERT(bed_input); Input* melting_input=inputs->GetInput(MeltingRateEnum); ISSMASSERT(melting_input); Input* vx_input=inputs->GetInput(VxEnum); ISSMASSERT(vx_input); @@ -3686,7 +3636,5 @@ vzstokes_input=inputs->GetInput(VzStokesEnum); ISSMASSERT(vzstokes_input); } - Input* bed_input=inputs->GetInput(BedEnum); ISSMASSERT(bed_input); - - /*For icesheets: */ + /* Start looping on the number of gaussian points: */ gauss=new GaussTria(2); @@ -3695,8 +3643,6 @@ gauss->GaussPoint(ig); - /*Get melting at gaussian point: */ melting_input->GetParameterValue(&meltingvalue, gauss); - - /*Get velocity at gaussian point: */ + bed_input->GetParameterDerivativeValue(&slope[0],&xyz_list[0][0],gauss); vx_input->GetParameterValue(&vx, gauss); vy_input->GetParameterValue(&vy, gauss); @@ -3706,26 +3652,17 @@ else vz=0; - /*Get bed slope: */ - bed_input->GetParameterDerivativeValue(&slope[0],&xyz_list[0][0],gauss); dbdx=slope[0]; dbdy=slope[1]; - /* Get Jacobian determinant: */ GetJacobianDeterminant3d(&Jdet, &xyz_list[0][0],gauss); - - //Get L matrix if viscous basal drag present: GetL(&L[0], &xyz_list[0][0], gauss,NDOF1); - /*Build gaussian vector: */ for(i=0;iweight*(vx*dbdx+vy*dbdy-vz-meltingvalue)*L[i]; } - /*Add pe_g_gaussian vector to pe_g: */ for( i=0; i