Index: /issm/trunk/src/c/objects/Loads/Icefront.cpp
===================================================================
--- /issm/trunk/src/c/objects/Loads/Icefront.cpp	(revision 5719)
+++ /issm/trunk/src/c/objects/Loads/Icefront.cpp	(revision 5720)
@@ -504,111 +504,27 @@
 void Icefront::CreatePVectorDiagnosticMacAyeal3d( Vec pg){
 
-	int       i;
-	const int numnodes= NUMVERTICESSEG;
-	const int numdofs = numnodes *NDOF2;
-	int*      doflist=NULL;
-	double    xyz_list[numnodes][3];
-
-	/*Objects: */
-	double    pe_g[numdofs] = {0.0};
-
-	/*Segment*/
-	double   normal[2];
-	int      num_gauss;
-	double*  segment_gauss_coord=NULL;
-	double*  gauss_weights=NULL;
-	double   gauss_weight;
-	double   gauss_coord;
-	int      ig;
-	double   Jdet;
-	double   L[2];
-	double   thickness,bed,pressure;
-	double   ice_pressure,water_pressure,air_pressure,rho_water,rho_ice,gravity;
-	double   surface_under_water,base_under_water;
-	int      fill;
-
-	/*Inputs*/
-	Input*  thickness_input=NULL;
-	Input*  bed_input=NULL;
-	Penta*  penta=NULL;
-
-	BeamRef* beam=NULL;
-
-	//check that the element is onbed (collapsed formulation) otherwise:pe=0
-	if(!element->GetOnBed()) return;
-
-	/*Recover inputs: */
-	thickness_input=((Penta*)element)->inputs->GetInput(ThicknessEnum);
-	bed_input      =((Penta*)element)->inputs->GetInput(BedEnum);
-	inputs->GetParameterValue(&fill,FillEnum);
-
-	/* Get dof list and node coordinates: */
-	GetDofList(&doflist,MacAyealApproximationEnum);
-	GetVerticesCoordinates(&xyz_list[0][0],nodes,numnodes);
-
-	/*Get Matpar parameters*/
-	rho_water=matpar->GetRhoWater();
-	rho_ice  =matpar->GetRhoIce();
-	gravity  =matpar->GetG();
-
-	/*Get normal*/
-	GetNormal(&normal[0],xyz_list);
-
-	//Get gaussian points and weights. order 2 since we have a product of 2 nodal functions
-	num_gauss=3;
-	gaussSegment(&segment_gauss_coord, &gauss_weights, num_gauss);
-
-	for(ig=0;ig<num_gauss;ig++){
-
-		/*Pick up the gaussian point: */
-		gauss_weight=gauss_weights[ig];
-		gauss_coord=segment_gauss_coord[ig];
-
-		//compute Jacobian of segment
-		beam->GetJacobianDeterminant2d(&Jdet,&xyz_list[0][0],gauss_coord);
-
-		/*Get thickness and bed*/
-		this->GetParameterValue(&thickness,gauss_coord,thickness_input);
-		this->GetParameterValue(&bed,      gauss_coord,bed_input);
-
-		/*Compute total pressure applied to ice front*/
-		if (fill==WaterEnum){
-			//icefront ends in water: 
-			ice_pressure=1.0/2.0*gravity*rho_ice*pow(thickness,2);
-			air_pressure=0;
-
-			//Now deal with water pressure
-			surface_under_water=min(0,thickness+bed); // 0 if the top of the glacier is above water level
-			base_under_water=min(0,bed);              // 0 if the bottom of the glacier is above water level
-			water_pressure=1.0/2.0*gravity*rho_water*(pow(surface_under_water,2) - pow(base_under_water,2));
-		}
-		else if (fill==AirEnum){
-			ice_pressure=1.0/2.0*gravity*rho_ice*pow(thickness,2);
-			air_pressure=0;
-			water_pressure=0;
-		}
-		else{
-			ISSMERROR("fill type %i not supported yet",fill);
-		}
-		pressure = ice_pressure + water_pressure + air_pressure;
-
-		/*Get L*/
-		beam->GetNodalFunctions(&L[0],gauss_coord);
-
-		for (i=0;i<numnodes;i++){
-			pe_g[2*i+0]+= pressure*Jdet*gauss_weights[ig]*normal[0]*L[i];
-			pe_g[2*i+1]+= pressure*Jdet*gauss_weights[ig]*normal[1]*L[i];
-		}
-	} //for(ig=0;ig<num_gauss;ig++)
-
-
-	/*Plug pe_g into vector: */
-	VecSetValues(pg,numdofs,doflist,pe_g,ADD_VALUES);
-
-	/*Free ressources:*/
-	xfree((void**)&segment_gauss_coord);
-	xfree((void**)&gauss_weights);
-	xfree((void**)&doflist);
-
+	Icefront* icefront=NULL;
+	Penta*    penta=NULL;
+	Tria*     tria=NULL;
+	bool      onbed;
+
+	/*Cast element onto Penta*/
+	penta   =(Penta*)this->element;
+
+	/*Return if not on bed*/
+	penta->inputs->GetParameterValue(&onbed,ElementOnBedEnum);
+	if(!onbed) return;
+
+	/*Spawn Tria and call MacAyeal2d*/
+	tria    =(Tria*)penta->SpawnTria(0,1,2);
+	icefront=(Icefront*)this->copy();
+	icefront->element=tria;
+	icefront->inputs->AddInput(new IntInput(TypeEnum,MacAyeal2dIceFrontEnum));
+	icefront->CreatePVectorDiagnosticMacAyeal2d(pg);
+
+	/*clean-up*/
+	delete tria->matice;
+	delete tria;
+	delete icefront;
 }
 /*}}}*/