Index: ../trunk-jpl/src/c/analyses/DamageEvolutionAnalysis.cpp
===================================================================
--- ../trunk-jpl/src/c/analyses/DamageEvolutionAnalysis.cpp	(revision 26222)
+++ ../trunk-jpl/src/c/analyses/DamageEvolutionAnalysis.cpp	(revision 26223)
@@ -808,6 +808,11 @@
 	delete gauss;
 	return Me;
 }/*}}}*/
+ElementVector* DamageEvolutionAnalysis::CreateFctPVector(Element* element){/*{{{*/
+
+	return this->CreatePVector(element);
+
+}/*}}}*/
 void           DamageEvolutionAnalysis::FctKMatrix(Matrix<IssmDouble>** pKff,Matrix<IssmDouble>** pKfs,FemModel* femmodel){/*{{{*/
 
 	/*Output*/
@@ -836,6 +841,26 @@
 		delete Kfs;
 	}
 }/*}}}*/
+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){/*{{{*/
 
 	/*Initialize Mass matrix*/
Index: ../trunk-jpl/src/c/analyses/DamageEvolutionAnalysis.h
===================================================================
--- ../trunk-jpl/src/c/analyses/DamageEvolutionAnalysis.h	(revision 26222)
+++ ../trunk-jpl/src/c/analyses/DamageEvolutionAnalysis.h	(revision 26223)
@@ -38,7 +38,9 @@
 		/*FCT*/
 		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: ../trunk-jpl/src/c/solutionsequences/solutionsequence_fct.cpp
===================================================================
--- ../trunk-jpl/src/c/solutionsequences/solutionsequence_fct.cpp	(revision 26222)
+++ ../trunk-jpl/src/c/solutionsequences/solutionsequence_fct.cpp	(revision 26223)
@@ -417,6 +417,7 @@
 			danalysis = new DamageEvolutionAnalysis();
 			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");
 	}
Index: ../trunk-jpl/src/c/classes/Elements/Tria.cpp
===================================================================
--- ../trunk-jpl/src/c/classes/Elements/Tria.cpp	(revision 26222)
+++ ../trunk-jpl/src/c/classes/Elements/Tria.cpp	(revision 26223)
@@ -1815,7 +1815,7 @@
 		total_weight+=gauss->weight;
 	}
 	//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;
 
 	/*free ressources:*/
Index: ../trunk-jpl/src/c/classes/SealevelGeometry.cpp
===================================================================
--- ../trunk-jpl/src/c/classes/SealevelGeometry.cpp	(revision 26222)
+++ ../trunk-jpl/src/c/classes/SealevelGeometry.cpp	(revision 26223)
@@ -19,7 +19,7 @@
 SealevelGeometry::SealevelGeometry(int localnelin){ /*{{{*/
 	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);
 		}
 		vlatbarycentre[i]=NULL; //we don't know yet 
@@ -42,7 +42,7 @@
 }; /*}}}*/
 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]);
 		}
 		xDelete<IssmDouble>(LoadArea[i]);
Index: ../trunk-jpl/src/c/classes/SealevelGeometry.h
===================================================================
--- ../trunk-jpl/src/c/classes/SealevelGeometry.h	(revision 26222)
+++ ../trunk-jpl/src/c/classes/SealevelGeometry.h	(revision 26223)
@@ -10,7 +10,7 @@
 #define SLGEOM_OCEAN 0 
 #define SLGEOM_ICE 1 
 #define SLGEOM_WATER 2
-#define MAXVERTICES 3
+#define SLMAXVERTICES 3
 
 #include "../toolkits/toolkits.h"
 
@@ -19,7 +19,7 @@
 	public: 
 
 		int         localnel;
-		IssmDouble* LoadWeigths[SLGEOM_NUMLOADS][MAXVERTICES];
+		IssmDouble* LoadWeigths[SLGEOM_NUMLOADS][SLMAXVERTICES];
 		IssmDouble* LoadArea[SLGEOM_NUMLOADS];
 		Vector<IssmDouble>* vlatbarycentre[SLGEOM_NUMLOADS];
 		Vector<IssmDouble>* vlongbarycentre[SLGEOM_NUMLOADS];