Index: /issm/trunk-jpl/src/c/analyses/DamageEvolutionAnalysis.cpp
===================================================================
--- /issm/trunk-jpl/src/c/analyses/DamageEvolutionAnalysis.cpp	(revision 26222)
+++ /issm/trunk-jpl/src/c/analyses/DamageEvolutionAnalysis.cpp	(revision 26223)
@@ -809,4 +809,9 @@
 	return Me;
 }/*}}}*/
+ElementVector* DamageEvolutionAnalysis::CreateFctPVector(Element* element){/*{{{*/
+
+	return this->CreatePVector(element);
+
+}/*}}}*/
 void           DamageEvolutionAnalysis::FctKMatrix(Matrix<IssmDouble>** pKff,Matrix<IssmDouble>** pKfs,FemModel* femmodel){/*{{{*/
 
@@ -837,4 +842,24 @@
 	}
 }/*}}}*/
+void           DamageEvolutionAnalysis::FctPVector(Vector<IssmDouble>** ppf,FemModel* femmodel){/*{{{*/
+
+	/*Output*/
+	Vector<IssmDouble>* pf = NULL;
+
+	/*Initialize P vector*/
+	AllocateSystemMatricesx(NULL,NULL,NULL,&pf,femmodel);
+
+	/*Create and assemble matrix*/
+	for(Object* & object : femmodel->elements->objects){
+		Element*       element = xDynamicCast<Element*>(object);
+		ElementVector* pe      = this->CreateFctPVector(element);
+		if(pe) pe->AddToGlobal(pf);
+		delete pe;
+	}
+	pf->Assemble();
+
+	/*Assign output pointer*/
+	*ppf=pf;
+}/*}}}*/
 void           DamageEvolutionAnalysis::MassMatrix(Matrix<IssmDouble>** pMff,FemModel* femmodel){/*{{{*/
 
Index: /issm/trunk-jpl/src/c/analyses/DamageEvolutionAnalysis.h
===================================================================
--- /issm/trunk-jpl/src/c/analyses/DamageEvolutionAnalysis.h	(revision 26222)
+++ /issm/trunk-jpl/src/c/analyses/DamageEvolutionAnalysis.h	(revision 26223)
@@ -39,6 +39,8 @@
 		ElementMatrix* CreateFctKMatrix(Element* element);
 		ElementMatrix* CreateMassMatrix(Element* element);
+		ElementVector* CreateFctPVector(Element* element);
 		void           FctKMatrix(Matrix<IssmDouble>** pKff,Matrix<IssmDouble>** pKfs,FemModel* femmodel);
 		void           MassMatrix(Matrix<IssmDouble>** pMff,FemModel* femmodel);
+		void				FctPVector(Vector<IssmDouble>** ppf,FemModel* femmodel);
 };
 #endif
Index: /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp	(revision 26222)
+++ /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp	(revision 26223)
@@ -1816,5 +1816,5 @@
 	}
 	//normalize to phi. 
-	if(total_weight)for(int i=0;i<NUMVERTICES;i++)weights[i]/=total_weight/phi; 
+	if(total_weight>0.) for(int i=0;i<NUMVERTICES;i++)weights[i]/=total_weight/phi; 
 	else for(int i=0;i<NUMVERTICES;i++)weights[i]=0;
 
Index: /issm/trunk-jpl/src/c/classes/SealevelGeometry.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/SealevelGeometry.cpp	(revision 26222)
+++ /issm/trunk-jpl/src/c/classes/SealevelGeometry.cpp	(revision 26223)
@@ -20,5 +20,5 @@
 	localnel=localnelin;
 	for(int i=0;i<SLGEOM_NUMLOADS;i++){
-		for (int j=0;j<MAXVERTICES;j++){
+		for (int j=0;j<SLMAXVERTICES;j++){
 			LoadWeigths[i][j]=xNewZeroInit<IssmDouble>(localnel);
 		}
@@ -43,5 +43,5 @@
 SealevelGeometry::~SealevelGeometry(){ /*{{{*/
 	for(int i=0;i<SLGEOM_NUMLOADS;i++){
-		for (int j=0;j<MAXVERTICES;j++){
+		for (int j=0;j<SLMAXVERTICES;j++){
 			xDelete<IssmDouble>(LoadWeigths[i][j]);
 		}
Index: /issm/trunk-jpl/src/c/classes/SealevelGeometry.h
===================================================================
--- /issm/trunk-jpl/src/c/classes/SealevelGeometry.h	(revision 26222)
+++ /issm/trunk-jpl/src/c/classes/SealevelGeometry.h	(revision 26223)
@@ -11,5 +11,5 @@
 #define SLGEOM_ICE 1 
 #define SLGEOM_WATER 2
-#define MAXVERTICES 3
+#define SLMAXVERTICES 3
 
 #include "../toolkits/toolkits.h"
@@ -20,5 +20,5 @@
 
 		int         localnel;
-		IssmDouble* LoadWeigths[SLGEOM_NUMLOADS][MAXVERTICES];
+		IssmDouble* LoadWeigths[SLGEOM_NUMLOADS][SLMAXVERTICES];
 		IssmDouble* LoadArea[SLGEOM_NUMLOADS];
 		Vector<IssmDouble>* vlatbarycentre[SLGEOM_NUMLOADS];
Index: /issm/trunk-jpl/src/c/solutionsequences/solutionsequence_fct.cpp
===================================================================
--- /issm/trunk-jpl/src/c/solutionsequences/solutionsequence_fct.cpp	(revision 26222)
+++ /issm/trunk-jpl/src/c/solutionsequences/solutionsequence_fct.cpp	(revision 26223)
@@ -418,4 +418,5 @@
 			danalysis->MassMatrix(&Mc,femmodel);
 			danalysis->FctKMatrix(&K,NULL,femmodel);
+			danalysis->FctPVector(&p,femmodel);
 			break;
 		default: _error_("analysis type " << EnumToStringx(analysis_type) << " not supported for FCT\n");