Index: /issm/trunk/src/c/EnumDefinitions/EnumDefinitions.h
===================================================================
--- /issm/trunk/src/c/EnumDefinitions/EnumDefinitions.h	(revision 4909)
+++ /issm/trunk/src/c/EnumDefinitions/EnumDefinitions.h	(revision 4910)
@@ -215,4 +215,5 @@
 	ResetPenaltiesEnum,
 	RheologyBEnum,
+	RheologyB2dEnum,
 	RheologyNEnum,
 	SegmentOnIceShelfEnum,
Index: /issm/trunk/src/c/Makefile.am
===================================================================
--- /issm/trunk/src/c/Makefile.am	(revision 4909)
+++ /issm/trunk/src/c/Makefile.am	(revision 4910)
@@ -153,4 +153,6 @@
 					./objects/Elements/Beam.h\
 					./objects/Elements/Beam.cpp\
+					./objects/Elements/BeamRef.h\
+					./objects/Elements/BeamRef.cpp\
 					./objects/Elements/Penta.h\
 					./objects/Elements/Penta.cpp\
@@ -709,4 +711,6 @@
 					./objects/Elements/Beam.h\
 					./objects/Elements/Beam.cpp\
+					./objects/Elements/BeamRef.h\
+					./objects/Elements/BeamRef.cpp\
 					./objects/Elements/Penta.h\
 					./objects/Elements/Penta.cpp\
Index: /issm/trunk/src/c/objects/Elements/Beam.cpp
===================================================================
--- /issm/trunk/src/c/objects/Elements/Beam.cpp	(revision 4909)
+++ /issm/trunk/src/c/objects/Elements/Beam.cpp	(revision 4910)
@@ -1058,31 +1058,4 @@
 }
 /*}}}*/
-/*FUNCTION Beam::GetJacobianDeterminant{{{1*/
-void Beam::GetJacobianDeterminant(double* pJdet,double* z_list, double gauss_coord){
-
-
-	double Jdet;
-
-	Jdet=1.0/2.0*(z_list[1]-z_list[0]);
-
-	if(Jdet<0){
-		ISSMERROR(" negative jacobian determinant!");
-	}
-	
-	*pJdet=Jdet;
-}
-/*}}}*/
-/*FUNCTION Beam::GetNodalFunctions{{{1*/
-void Beam::GetNodalFunctions(double* l1l2, double gauss_coord){
-	
-	/*This routine returns the values of the nodal functions  at the gaussian point.*/
-
-	/*First nodal function: */
-	l1l2[0]=-0.5*gauss_coord+0.5;
-
-	/*Second nodal function: */
-	l1l2[1]=+0.5*gauss_coord+0.5;
-}
-/*}}}*/
 /*FUNCTION Beam::GetParameterValue{{{1*/
 void Beam::GetParameterValue(double* pvalue, double* value_list,double gauss_coord){
Index: /issm/trunk/src/c/objects/Elements/Beam.h
===================================================================
--- /issm/trunk/src/c/objects/Elements/Beam.h	(revision 4909)
+++ /issm/trunk/src/c/objects/Elements/Beam.h	(revision 4910)
@@ -10,4 +10,5 @@
 
 #include "./Element.h"
+#include "./BeamRef.h"
 class Hook;
 class Parameters;
@@ -22,5 +23,5 @@
 /*}}}*/
 
-class Beam: public Element{
+class Beam: public Element,public BeamRef {
 
 	public:
@@ -67,8 +68,8 @@
 		void	   Configure(Elements* elements,Loads* loads,DataSet* nodes,Materials* materials,Parameters* parameters);
 		void	   SetCurrentConfiguration(Elements* elements,Loads* loads,DataSet* nodes,Materials* materials,Parameters* parameters);
-		double	   CostFunction(void);
+		double	CostFunction(void);
 		void	   CreateKMatrix(Mat Kgg);
 		void	   CreatePVector(Vec pg);
-		void   DeleteResults(void);
+		void     DeleteResults(void);
 		void	   GetBedList(double* bed_list);
 		void*	   GetMatPar();
@@ -102,11 +103,11 @@
 		void	   MinVy(double* pminvy, bool process_units);
 		void	   MinVz(double* pminvz, bool process_units);
-		double	   Misfit(void);
+		double	Misfit(void);
 		void	   PatchFill(int* pcount, Patch* patch);
 		void	   PatchSize(int* pnumrows, int* pnumvertices,int* pnumnodes);
 		void	   ProcessResultsUnits(void);
-		double	   SurfaceArea(void);
-		void       Update(int index, IoModel* iomodel,int analysis_counter,int analysis_type);
-		void       UpdateGeometry(void){ISSMERROR("not implemented yet");};
+		double	SurfaceArea(void);
+		void     Update(int index, IoModel* iomodel,int analysis_counter,int analysis_type);
+		void     UpdateGeometry(void){ISSMERROR("not implemented yet");};
 		/*}}}*/
 		/*Beam specific routines: {{{1*/
@@ -115,6 +116,4 @@
 		void	  GetDofList(int* doflist,int* pnumberofdofs);
 		void	  GetDofList1(int* doflist);
-		void	  GetJacobianDeterminant(double* pJdet,double* z_list, double gauss_coord);
-		void	  GetNodalFunctions(double* l1l2, double gauss_coord);
 		void	  GetParameterValue(double* pp, double* plist, double* gauss_l1l2l3);
 		void	  GetParameterValue(double* pvalue, double* value_list,double gauss_coord);
Index: /issm/trunk/src/c/objects/Elements/BeamRef.cpp
===================================================================
--- /issm/trunk/src/c/objects/Elements/BeamRef.cpp	(revision 4910)
+++ /issm/trunk/src/c/objects/Elements/BeamRef.cpp	(revision 4910)
@@ -0,0 +1,70 @@
+/*!\file BeamRef.c
+ * \brief: implementation of the BeamRef object
+ */
+
+/*Headers:*/
+/*{{{1*/
+#ifdef HAVE_CONFIG_H
+	#include "config.h"
+#else
+#error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!"
+#endif
+
+#include "stdio.h"
+#include <string.h>
+#include "../objects.h"
+#include "../../EnumDefinitions/EnumDefinitions.h"
+#include "../../shared/shared.h"
+#include "../../Container/Container.h"
+#include "../../include/include.h"
+/*}}}*/
+
+/*Object constructors and destructor*/
+/*FUNCTION BeamRef::BeamRef(){{{1*/
+BeamRef::BeamRef(){
+	this->element_type_list=NULL;
+}
+/*}}}*/
+/*FUNCTION BeamRef::BeamRef(int* types,int nummodels){{{1*/
+BeamRef::BeamRef(const int nummodels){
+
+	/*Only allocate pointer*/
+	element_type_list=(int*)xmalloc(nummodels*sizeof(int));
+
+}
+/*}}}*/
+/*FUNCTION BeamRef::~BeamRef(){{{1*/
+BeamRef::~BeamRef(){
+	xfree((void**)&element_type_list);
+}
+/*}}}*/
+
+/*Management*/
+/*FUNCTION BeamRef::SetElementType{{{1*/
+void BeamRef::SetElementType(int type,int type_counter){
+
+	ISSMASSERT(type==P1Enum || type==P1DGEnum);
+
+	/*initialize element type*/
+	this->element_type_list[type_counter]=type;
+}
+/*}}}*/
+
+/*Reference Element numerics*/
+/*FUNCTION BeamRef::GetJacobianDeterminant {{{1*/
+void BeamRef::GetJacobianDeterminant(double*  Jdet, double* z_list,double gauss){
+
+	*Jdet=1.0/2.0*(z_list[1]-z_list[0]);
+	if(*Jdet<0) ISSMERROR(" negative jacobian determinant!");
+
+}
+/*}}}*/
+/*FUNCTION BeamRef::GetNodalFunctions {{{1*/
+void BeamRef::GetNodalFunctions(double* l1l2, double gauss){
+	/*This routine returns the values of the nodal functions  at the gaussian point.*/
+
+	l1l2[0]=-0.5*gauss+0.5;
+	l1l2[1]=+0.5*gauss+0.5;
+
+}
+/*}}}*/
Index: /issm/trunk/src/c/objects/Elements/BeamRef.h
===================================================================
--- /issm/trunk/src/c/objects/Elements/BeamRef.h	(revision 4910)
+++ /issm/trunk/src/c/objects/Elements/BeamRef.h	(revision 4910)
@@ -0,0 +1,29 @@
+/*!\file:  BeamRef.h
+ * \brief abstract class for handling Beam oriented routines, like nodal functions, 
+ * strain rate generation, etc ...
+ */ 
+
+
+#ifndef _BEAMREF_H_
+#define _BEAMREF_H_
+
+class BeamRef{
+	
+
+	public: 
+		int* element_type_list; //P1CG, P1DG, MINI, P2...
+		int  element_type;
+		
+		BeamRef();
+		BeamRef(const int nummodels);
+		~BeamRef();
+
+		/*Management*/
+		void SetElementType(int type,int type_counter);
+
+		/*Numerics*/
+		void GetJacobianDeterminant(double* Jdet, double* z_list,double gauss);
+		void GetNodalFunctions(double* l1l2, double gauss);
+
+};
+#endif
Index: /issm/trunk/src/c/objects/Elements/TriaRef.cpp
===================================================================
--- /issm/trunk/src/c/objects/Elements/TriaRef.cpp	(revision 4909)
+++ /issm/trunk/src/c/objects/Elements/TriaRef.cpp	(revision 4910)
@@ -256,5 +256,5 @@
 	/*Get Determinant*/
 	Matrix2x2Determinant(Jdet,&J[0][0]);
-	if(Jdet<0) ISSMERROR("negative jacobian determinant!");
+	if(*Jdet<0) ISSMERROR("negative jacobian determinant!");
 
 }
@@ -277,10 +277,6 @@
 	z3=*(xyz_list+3*2+2);
 
-
 	*Jdet=SQRT3/6.0*pow(pow(((y2-y1)*(z3-z1)-(z2-z1)*(y3-y1)),2.0)+pow(((z2-z1)*(x3-x1)-(x2-x1)*(z3-z1)),2.0)+pow(((x2-x1)*(y3-y1)-(y2-y1)*(x3-x1)),2.0),0.5);
-
-	if(Jdet<0){
-		ISSMERROR("negative jacobian determinant!");
-	}
+	if(*Jdet<0) ISSMERROR("negative jacobian determinant!");
 
 }
Index: /issm/trunk/src/c/objects/Inputs/BeamVertexInput.h
===================================================================
--- /issm/trunk/src/c/objects/Inputs/BeamVertexInput.h	(revision 4909)
+++ /issm/trunk/src/c/objects/Inputs/BeamVertexInput.h	(revision 4910)
@@ -10,8 +10,9 @@
 /*{{{1*/
 #include "./Input.h"
+#include "../Elements/BeamRef.h"
 /*}}}*/
 
 
-class BeamVertexInput: public Input{
+class BeamVertexInput: public Input,public BeamRef{
 
 	public:
Index: /issm/trunk/src/c/objects/objects.h
===================================================================
--- /issm/trunk/src/c/objects/objects.h	(revision 4909)
+++ /issm/trunk/src/c/objects/objects.h	(revision 4910)
@@ -34,4 +34,5 @@
 /*Elements: */
 #include "./Elements/Beam.h"
+#include "./Elements/BeamRef.h"
 #include "./Elements/Element.h"
 #include "./Elements/Penta.h"