Index: /issm/trunk/src/c/objects/Tria.cpp
===================================================================
--- /issm/trunk/src/c/objects/Tria.cpp	(revision 599)
+++ /issm/trunk/src/c/objects/Tria.cpp	(revision 600)
@@ -579,5 +579,5 @@
 
 	/*recover extra inputs from users, at current convergence iteration: */
-	found=inputs->Recover("velocity",&vxvy_list[0][0],2,dofs,numgrids,(void**)nodes);
+	found=inputs->Recover("velocity_average",&vxvy_list[0][0],2,dofs,numgrids,(void**)nodes);
 	if(!found)throw ErrorException(__FUNCT__," could not find velocity_average  in inputs!");
 
@@ -728,4 +728,102 @@
 }
 
+#undef __FUNCT__ 
+#define __FUNCT__ "Tria::CreatePVectorPrognostic"
+void  Tria::CreatePVectorPrognostic(Vec pg ,void* vinputs,int analysis_type,int sub_analysis_type){
+
+
+	/* local declarations */
+	int             i,j;
+
+	/* node data: */
+	const int    numgrids=3;
+	const int    NDOF1=1;
+	const int    numdof=NDOF1*numgrids;
+	double       xyz_list[numgrids][3];
+	int          doflist[numdof];
+	int          numberofdofspernode;
+	
+	/* gaussian points: */
+	int     num_gauss,ig;
+	double* first_gauss_area_coord  =  NULL;
+	double* second_gauss_area_coord =  NULL;
+	double* third_gauss_area_coord  =  NULL;
+	double* gauss_weights           =  NULL;
+	double  gauss_weight;
+	double  gauss_l1l2l3[3];
+
+	/* matrix */
+	double pe_g[numgrids]={0.0};
+	double L[numgrids];
+	double Jdettria;
+	
+	/*input parameters for structural analysis (diagnostic): */
+	double  accumulation_list[numgrids]={0.0};
+	double  accumulation_g;
+	double  melting_list[numgrids]={0.0};
+	double  melting_g;
+	double  thickness_list[numgrids]={0.0};
+	double  thickness_g;
+	double  dt;
+	int     dofs[1]={1};
+	int     found=0;
+
+	ParameterInputs* inputs=NULL;
+
+	/*recover pointers: */
+	inputs=(ParameterInputs*)vinputs;
+
+	/*recover extra inputs from users, at current convergence iteration: */
+	found=inputs->Recover("accumulation",&accumulation_list[0],1,dofs,numgrids,(void**)nodes);
+	if(!found)throw ErrorException(__FUNCT__," could not find accumulation in inputs!");
+	found=inputs->Recover("melting",&melting_list[0],1,dofs,numgrids,(void**)nodes);
+	if(!found)throw ErrorException(__FUNCT__," could not find melting in inputs!");
+	found=inputs->Recover("thickness",&thickness_list[0],1,dofs,numgrids,(void**)nodes);
+	if(!found)throw ErrorException(__FUNCT__," could not find thickness in inputs!");
+	found=inputs->Recover("dt",&dt);
+	if(!found)throw ErrorException(__FUNCT__," could not find dt in inputs!");
+
+	/* Get node coordinates and dof list: */
+	GetElementNodeData( &xyz_list[0][0], nodes, numgrids);
+	GetDofList(&doflist[0],&numberofdofspernode);
+
+	/* Get gaussian points and weights (make this a statically initialized list of points? fstd): */
+	GaussTria( &num_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &gauss_weights, 2);
+
+	/* Start  looping on the number of gaussian points: */
+	for (ig=0; ig<num_gauss; ig++){
+		/*Pick up the gaussian point: */
+		gauss_weight=*(gauss_weights+ig);
+		gauss_l1l2l3[0]=*(first_gauss_area_coord+ig); 
+		gauss_l1l2l3[1]=*(second_gauss_area_coord+ig);
+		gauss_l1l2l3[2]=*(third_gauss_area_coord+ig);
+
+		/* Get Jacobian determinant: */
+		GetJacobianDeterminant2d(&Jdettria, &xyz_list[0][0],gauss_l1l2l3);
+
+		/*Get L matrix: */
+		GetL(&L[0], &xyz_list[0][0], gauss_l1l2l3,numberofdofspernode);
+
+		/* Get accumulation, melting and thickness at gauss point */
+		GetParameterValue(&accumulation_g, &accumulation_list[0],gauss_l1l2l3);
+		GetParameterValue(&melting_g, &melting_list[0],gauss_l1l2l3);
+		GetParameterValue(&thickness_g, &thickness_list[0],gauss_l1l2l3);
+		
+		/* Add value into pe_g: */
+		for( i=0; i<numdof; i++) pe_g[i]+=Jdettria*gauss_weight*(thickness_g+dt*(accumulation_g-melting_g))*L[i];
+		
+	} // for (ig=0; ig<num_gauss; ig++)
+
+	/*Add pe_g to global matrix Kgg: */
+	VecSetValues(pg,numdof,doflist,(const double*)pe_g,ADD_VALUES);
+
+	cleanup_and_return: 
+	xfree((void**)&first_gauss_area_coord);
+	xfree((void**)&second_gauss_area_coord);
+	xfree((void**)&third_gauss_area_coord);
+	xfree((void**)&gauss_weights);
+
+}
+
 
 #undef __FUNCT__ 
@@ -1000,4 +1098,8 @@
 		
 		CreatePVectorSlopeCompute( pg,inputs,analysis_type,sub_analysis_type);
+	}
+	else if (analysis_type==PrognosticAnalysisEnum()){
+
+		CreatePVectorPrognostic( pg,inputs,analysis_type,sub_analysis_type);
 	}
 	else{
Index: /issm/trunk/src/c/objects/Tria.h
===================================================================
--- /issm/trunk/src/c/objects/Tria.h	(revision 599)
+++ /issm/trunk/src/c/objects/Tria.h	(revision 600)
@@ -115,4 +115,5 @@
 		void  CreatePVectorThermalSheet( Vec pg, void* vinputs, int analysis_type,int sub_analysis_type);
 		void  CreateKMatrixPrognostic(Mat Kgg,void* vinputs,int analysis_type,int sub_analysis_type);
+		void  CreatePVectorPrognostic(Vec pg,void* vinputs,int analysis_type,int sub_analysis_type);