Index: /issm/trunk/src/c/objects/Beam.cpp
===================================================================
--- /issm/trunk/src/c/objects/Beam.cpp	(revision 2906)
+++ /issm/trunk/src/c/objects/Beam.cpp	(revision 2907)
@@ -18,12 +18,16 @@
 
 
+/*Object constructors and destructor*/
+/*FUNCTION Beam constructor {{{1*/
 Beam::Beam(){
 	return;
 }
-
+/*}}}*/
+/*FUNCTION Beam destructor {{{1*/
 Beam::~Beam(){
 	return;
 }
-		
+/*}}}*/
+/*FUNCTION Beam creation {{{1*/
 Beam::Beam(int beam_id, int beam_mid, int beam_mparid, int beam_numparid, int beam_g[2], double beam_h[2], double beam_s[2],double beam_b[2],double beam_k[2],bool beam_onbed){
 
@@ -53,56 +57,8 @@
 	return;
 }
-
-#undef __FUNCT__
-#define __FUNCT__ "Beam::Echo"
-void Beam::Echo(void){
-
-	printf("Beam:\n");
-	printf("   id: %i\n",id);
-	printf("   mid: %i\n",mid);
-	printf("   mparid: %i\n",mparid);
-	printf("   node_ids=[%i,%i]\n",node_ids[0],node_ids[1]);
-	printf("   node_offsets=[%i,%i]\n",node_offsets[0],node_offsets[1]);
-	printf("   matice_offset=%i\n",matice_offset);
-	printf("   matpar_offset=%i\n",matpar_offset);
-	printf("   h=[%g,%g]\n",h[0],h[1]);
-	printf("   s=[%g,%g]\n",s[0],s[1]);
-	printf("   b=[%g,%g]\n",b[0],b[1]);
-	printf("   k=[%g,%g]\n",k[0],k[1]);
-	printf("   onbed=%i\n",onbed);
-	printf("   nodes: \n");
-	if(nodes[0])nodes[0]->Echo();
-	if(nodes[1])nodes[1]->Echo();
-	if(matice)matice->Echo();
-	if(matpar)matpar->Echo();
-
-	return;
-}
-
-#undef __FUNCT__
-#define __FUNCT__ "Beam::DeepEcho"
-void Beam::DeepEcho(void){
-
-	printf("Beam:\n");
-	printf("   id: %i\n",id);
-	printf("   mid: %i\n",mid);
-	printf("   mparid: %i\n",mparid);
-	printf("   node_ids=[%i,%i]\n",node_ids[0],node_ids[1]);
-	printf("   node_offsets=[%i,%i]\n",node_offsets[0],node_offsets[1]);
-	printf("   matice_offset=%i\n",matice_offset);
-	printf("   matpar_offset=%i\n",matpar_offset);
-	printf("   h=[%g,%g]\n",h[0],h[1]);
-	printf("   s=[%g,%g]\n",s[0],s[1]);
-	printf("   b=[%g,%g]\n",b[0],b[1]);
-	printf("   k=[%g,%g]\n",k[0],k[1]);
-	printf("   onbed=%i\n",onbed);
-	printf("   nodes: \n");
-	if(nodes[0])nodes[0]->Echo();
-	if(nodes[1])nodes[1]->Echo();
-	if(matice)matice->Echo();
-	if(matpar)matpar->Echo();
-
-	return;
-}
+/*}}}*/
+
+/*Object marshall*/
+/*FUNCTION Beam Marshall{{{1*/
 void  Beam::Marshall(char** pmarshalled_dataset){
 
@@ -142,5 +98,6 @@
 	return;
 }
-		
+/*}}}*/
+/*FUNCTION Beam MarshallSize{{{1*/
 int   Beam::MarshallSize(){
 	return sizeof(id)
@@ -164,9 +121,6 @@
 		+sizeof(int); //sizeof(int) for enum type
 }
-
-char* Beam::GetName(void){
-	return "beam";
-}
-
+/*}}}*/
+/*FUNCTION Beam Demarshall{{{1*/
 void  Beam::Demarshall(char** pmarshalled_dataset){
 
@@ -209,17 +163,38 @@
 	return;
 }
-int Beam::Enum(void){
-
-	return BeamEnum();
-
-}
-int    Beam::GetId(void){ return id; }
-
-int    Beam::MyRank(void){ 
-	extern int my_rank;
-	return my_rank; 
-}
-
-
+/*}}}*/
+
+/*Object functions*/
+/*FUNCTION Beam ComputePressure{{{1*/
+#undef __FUNCT__ 
+#define __FUNCT__ "Beam::ComputePressure"
+void  Beam::ComputePressure(Vec p_g){
+
+	int i;
+	const int numgrids=2;
+	int doflist[numgrids];
+	double pressure[numgrids];
+	double rho_ice,g;
+	double xyz_list[numgrids][3];
+
+	/*Get node data: */
+	GetElementNodeData( &xyz_list[0][0], nodes, numgrids);
+        
+	/*Get dof list on which we will plug the pressure values: */
+	GetDofList1(&doflist[0]);
+
+	/*pressure is lithostatic: */
+	rho_ice=matpar->GetRhoIce();
+	g=matpar->GetG();
+	for(i=0;i<numgrids;i++){
+		pressure[i]=rho_ice*g*(s[i]-xyz_list[i][2]);
+	}
+	
+	/*plug local pressure values into global pressure vector: */
+	VecSetValues(p_g,numgrids,doflist,(const double*)pressure,INSERT_VALUES);
+
+}
+/*}}}*/
+/*FUNCTION Beam Configure{{{1*/
 #undef __FUNCT__ 
 #define __FUNCT__ "Beam::Configure"
@@ -250,5 +225,13 @@
 
 }
-
+/*}}}*/
+/*FUNCTION Beam copy{{{1*/
+Object* Beam::copy() {
+	
+	return new Beam(*this); 
+
+}
+/*}}}*/
+/*FUNCTION Beam CreateKMatrix{{{1*/
 #undef __FUNCT__ 
 #define __FUNCT__ "Beam::CreateKMatrix"
@@ -271,6 +254,6 @@
 
 }
-
-
+/*}}}*/
+/*FUNCTION Beam CreateKMatrixDiagnosticHutter{{{1*/
 #undef __FUNCT__ 
 #define __FUNCT__ "Beam::CreateKMatrixDiagnosticHutter"
@@ -307,6 +290,6 @@
 
 }
-
-
+/*}}}*/
+/*FUNCTION Beam CreatePVector{{{1*/
 #undef __FUNCT__ 
 #define __FUNCT__ "Beam::CreatePVector"
@@ -326,6 +309,6 @@
 
 }
-
-
+/*}}}*/
+/*FUNCTION Beam CreatePVectorDiagnosticHutter{{{1*/
 #undef __FUNCT__ 
 #define __FUNCT__ "Beam::CreatePVectorDiagnosticHutter"
@@ -442,6 +425,270 @@
 	xfree((void**)&gauss_weights);
 }
-
-
+/*}}}*/
+/*FUNCTION Beam DeepEcho{{{1*/
+#undef __FUNCT__
+#define __FUNCT__ "Beam::DeepEcho"
+void Beam::DeepEcho(void){
+
+	printf("Beam:\n");
+	printf("   id: %i\n",id);
+	printf("   mid: %i\n",mid);
+	printf("   mparid: %i\n",mparid);
+	printf("   node_ids=[%i,%i]\n",node_ids[0],node_ids[1]);
+	printf("   node_offsets=[%i,%i]\n",node_offsets[0],node_offsets[1]);
+	printf("   matice_offset=%i\n",matice_offset);
+	printf("   matpar_offset=%i\n",matpar_offset);
+	printf("   h=[%g,%g]\n",h[0],h[1]);
+	printf("   s=[%g,%g]\n",s[0],s[1]);
+	printf("   b=[%g,%g]\n",b[0],b[1]);
+	printf("   k=[%g,%g]\n",k[0],k[1]);
+	printf("   onbed=%i\n",onbed);
+	printf("   nodes: \n");
+	if(nodes[0])nodes[0]->Echo();
+	if(nodes[1])nodes[1]->Echo();
+	if(matice)matice->Echo();
+	if(matpar)matpar->Echo();
+
+	return;
+}
+/*}}}*/
+/*FUNCTION Beam Du{{{1*/
+#undef __FUNCT__ 
+#define __FUNCT__ "Beam::Du"
+void  Beam::Du(_p_Vec*,void*,int,int){
+	throw ErrorException(__FUNCT__," not supported yet!");
+}
+/*}}}*/
+/*FUNCTION Beam Echo {{{1*/
+#undef __FUNCT__
+#define __FUNCT__ "Beam::Echo"
+void Beam::Echo(void){
+
+	printf("Beam:\n");
+	printf("   id: %i\n",id);
+	printf("   mid: %i\n",mid);
+	printf("   mparid: %i\n",mparid);
+	printf("   node_ids=[%i,%i]\n",node_ids[0],node_ids[1]);
+	printf("   node_offsets=[%i,%i]\n",node_offsets[0],node_offsets[1]);
+	printf("   matice_offset=%i\n",matice_offset);
+	printf("   matpar_offset=%i\n",matpar_offset);
+	printf("   h=[%g,%g]\n",h[0],h[1]);
+	printf("   s=[%g,%g]\n",s[0],s[1]);
+	printf("   b=[%g,%g]\n",b[0],b[1]);
+	printf("   k=[%g,%g]\n",k[0],k[1]);
+	printf("   onbed=%i\n",onbed);
+	printf("   nodes: \n");
+	if(nodes[0])nodes[0]->Echo();
+	if(nodes[1])nodes[1]->Echo();
+	if(matice)matice->Echo();
+	if(matpar)matpar->Echo();
+
+	return;
+}
+/*}}}*/
+/*FUNCTION Beam Enum{{{1*/
+int Beam::Enum(void){
+
+	return BeamEnum();
+
+}
+/*}}}*/
+/*FUNCTION Beam GetBedList{{{1*/
+#undef __FUNCT__ 
+#define __FUNCT__ "Beam::GetBedList"
+void  Beam::GetBedList(double*){
+	throw ErrorException(__FUNCT__," not supported yet!");
+}
+/*}}}*/
+/*FUNCTION Beam GetDofList{{{1*/
+void  Beam::GetDofList(int* doflist,int* pnumberofdofspernode){
+
+	int i,j;
+	int doflist_per_node[MAXDOFSPERNODE];
+	int numberofdofspernode;
+	
+	for(i=0;i<2;i++){
+		nodes[i]->GetDofList(&doflist_per_node[0],&numberofdofspernode);
+		for(j=0;j<numberofdofspernode;j++){
+			doflist[i*numberofdofspernode+j]=doflist_per_node[j];
+		}
+	}
+
+	/*Assign output pointers:*/
+	*pnumberofdofspernode=numberofdofspernode;
+
+}
+/*}}}*/
+/*FUNCTION Beam GetDofList1{{{1*/
+void  Beam::GetDofList1(int* doflist){
+
+	int i;
+	for(i=0;i<2;i++){
+		doflist[i]=nodes[i]->GetDofList1();
+	}
+
+}
+/*}}}*/
+/*FUNCTION Beam GetId{{{1*/
+int    Beam::GetId(void){ return id; }
+/*}}}*/
+/*FUNCTION Beam GetJacobianDeterminant{{{1*/
+#undef __FUNCT__ 
+#define __FUNCT__ "Beam::GetJacobianDeterminant"
+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){
+		throw ErrorException(__FUNCT__," negative jacobian determinant!");
+	}
+	
+	*pJdet=Jdet;
+}
+/*}}}*/
+/*FUNCTION Beam GetMatPar{{{1*/
+void* Beam::GetMatPar(){
+	return matpar;
+}
+/*}}}*/
+/*FUNCTION Beam GetName{{{1*/
+char* Beam::GetName(void){
+	return "beam";
+}
+/*}}}*/
+/*FUNCTION Beam GetNodalFunctions{{{1*/
+#undef __FUNCT__ 
+#define __FUNCT__ "Beam::GetNodalFunctions"
+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]=.5*gauss_coord+.5;
+
+	/*Second nodal function: */
+	l1l2[1]=-.5*gauss_coord+.5;
+}
+/*}}}*/
+/*FUNCTION Beam GetNodes{{{1*/
+void  Beam::GetNodes(void** vpnodes){
+	int i;
+	Node** pnodes=(Node**)vpnodes;
+
+	for(i=0;i<3;i++){
+		pnodes[i]=nodes[i];
+	}
+}
+/*}}}*/
+/*FUNCTION Beam GetOnBed{{{1*/
+#undef __FUNCT__ 
+#define __FUNCT__ "Beam::GetOnBed"
+int   Beam::GetOnBed(){
+	throw ErrorException(__FUNCT__," not supported yet!");
+}
+/*}}}*/
+/*FUNCTION Beam GetParameterValue{{{1*/
+#undef __FUNCT__ 
+#define __FUNCT__ "Beam::GetParameterValue"
+void Beam::GetParameterValue(double* pvalue, double* value_list,double gauss_coord){
+
+	double l1l2[2];
+	
+	GetNodalFunctions(&l1l2[0],gauss_coord);
+
+	*pvalue=l1l2[0]*value_list[0]+l1l2[1]*value_list[1];
+}
+/*}}}*/
+/*FUNCTION Beam GetShelf{{{1*/
+#undef __FUNCT__ 
+#define __FUNCT__ "Beam::GetShelf"
+int   Beam::GetShelf(){
+	throw ErrorException(__FUNCT__," not supported yet!");
+}
+/*}}}*/
+/*FUNCTION Beam GetThicknessList{{{1*/
+#undef __FUNCT__ 
+#define __FUNCT__ "Beam::GetThicknessList"
+void  Beam::GetThicknessList(double* thickness_list){
+	throw ErrorException(__FUNCT__," not supported yet!");
+}
+/*}}}*/
+/*FUNCTION Beam Gradj{{{1*/
+#undef __FUNCT__ 
+#define __FUNCT__ "Beam::Gradj"
+void  Beam::Gradj(_p_Vec*, void*, int, int,char*){
+	throw ErrorException(__FUNCT__," not supported yet!");
+}
+/*}}}*/
+/*FUNCTION Beam GradjB{{{1*/
+#undef __FUNCT__ 
+#define __FUNCT__ "Beam::GradjB"
+void  Beam::GradjB(_p_Vec*, void*, int, int){
+	throw ErrorException(__FUNCT__," not supported yet!");
+}
+/*}}}*/
+/*FUNCTION Beam GradjDrag{{{1*/
+#undef __FUNCT__ 
+#define __FUNCT__ "Beam::GradjDrag"
+void  Beam::GradjDrag(_p_Vec*, void*, int,int ){
+	throw ErrorException(__FUNCT__," not supported yet!");
+}
+/*}}}*/
+/*FUNCTION Beam MassFlux{{{1*/
+#undef __FUNCT__ 
+#define __FUNCT__ "Beam::MassFlux"
+double Beam::MassFlux( double* segment,double* ug){
+	throw ErrorException(__FUNCT__," not supported yet!");
+}
+/*}}}*/
+/*FUNCTION Beam MaticeConfiguration{{{1*/
+#undef __FUNCT__ 
+#define __FUNCT__ "Beam::MaticeConfiguration"
+void  Beam::MaticeConfiguration(Matice* beam_matice,int beam_matice_offset){
+	matice=beam_matice;
+	matice_offset=beam_matice_offset;
+}
+/*}}}*/
+/*FUNCTION Beam MatparConfiguration{{{1*/
+#undef __FUNCT__ 
+#define __FUNCT__ "Beam::MatparConfiguration"
+void  Beam::MatparConfiguration(Matpar* beam_matpar,int beam_matpar_offset){
+
+	matpar=beam_matpar;
+	matpar_offset=beam_matpar_offset;
+
+}
+/*}}}*/
+/*FUNCTION Beam Misfit{{{1*/
+#undef __FUNCT__ 
+#define __FUNCT__ "Beam::Misfit"
+double Beam::Misfit(void*,int,int ){
+	throw ErrorException(__FUNCT__," not supported yet!");
+}
+/*}}}*/
+/*FUNCTION Beam MyRank{{{1*/
+int    Beam::MyRank(void){ 
+	extern int my_rank;
+	return my_rank; 
+}
+/*}}}*/
+/*FUNCTION Beam NodeConfiguration{{{1*/
+#undef __FUNCT__ 
+#define __FUNCT__ "Beam::NodeConfiguration"
+void  Beam::NodeConfiguration(int* beam_node_ids,Node* beam_nodes[2],int* beam_node_offsets){
+
+	int i;
+	for(i=0;i<2;i++){
+		node_ids[i]=beam_node_ids[i];
+		nodes[i]=beam_nodes[i];
+		node_offsets[i]=beam_node_offsets[i];
+	}
+
+}
+/*}}}*/
+/*FUNCTION Beam UpdateFromInputs{{{1*/
 #undef __FUNCT__ 
 #define __FUNCT__ "Beam::UpdateFromInputs"
@@ -483,207 +730,3 @@
 
 }
-		
-void  Beam::GetDofList(int* doflist,int* pnumberofdofspernode){
-
-	int i,j;
-	int doflist_per_node[MAXDOFSPERNODE];
-	int numberofdofspernode;
-	
-	for(i=0;i<2;i++){
-		nodes[i]->GetDofList(&doflist_per_node[0],&numberofdofspernode);
-		for(j=0;j<numberofdofspernode;j++){
-			doflist[i*numberofdofspernode+j]=doflist_per_node[j];
-		}
-	}
-
-	/*Assign output pointers:*/
-	*pnumberofdofspernode=numberofdofspernode;
-
-}
-
-void  Beam::GetDofList1(int* doflist){
-
-	int i;
-	for(i=0;i<2;i++){
-		doflist[i]=nodes[i]->GetDofList1();
-	}
-
-}
-
-
-void* Beam::GetMatPar(){
-	return matpar;
-}
-
-		
-void  Beam::GetNodes(void** vpnodes){
-	int i;
-	Node** pnodes=(Node**)vpnodes;
-
-	for(i=0;i<3;i++){
-		pnodes[i]=nodes[i];
-	}
-}
-		
-
-Object* Beam::copy() {
-	
-	return new Beam(*this); 
-
-}
-
-
-#undef __FUNCT__ 
-#define __FUNCT__ "Beam::NodeConfiguration"
-void  Beam::NodeConfiguration(int* beam_node_ids,Node* beam_nodes[2],int* beam_node_offsets){
-
-	int i;
-	for(i=0;i<2;i++){
-		node_ids[i]=beam_node_ids[i];
-		nodes[i]=beam_nodes[i];
-		node_offsets[i]=beam_node_offsets[i];
-	}
-
-}
-#undef __FUNCT__ 
-#define __FUNCT__ "Beam::MaticeConfiguration"
-void  Beam::MaticeConfiguration(Matice* beam_matice,int beam_matice_offset){
-	matice=beam_matice;
-	matice_offset=beam_matice_offset;
-}
-
-#undef __FUNCT__ 
-#define __FUNCT__ "Beam::MatparConfiguration"
-void  Beam::MatparConfiguration(Matpar* beam_matpar,int beam_matpar_offset){
-
-	matpar=beam_matpar;
-	matpar_offset=beam_matpar_offset;
-
-}
-
-
-
-#undef __FUNCT__ 
-#define __FUNCT__ "Beam::GetShelf"
-int   Beam::GetShelf(){
-	throw ErrorException(__FUNCT__," not supported yet!");
-}
-#undef __FUNCT__ 
-#define __FUNCT__ "Beam::GetOnBed"
-int   Beam::GetOnBed(){
-	throw ErrorException(__FUNCT__," not supported yet!");
-}
-
-
-#undef __FUNCT__ 
-#define __FUNCT__ "Beam::GetBedList"
-void  Beam::GetBedList(double*){
-	throw ErrorException(__FUNCT__," not supported yet!");
-}
-
-#undef __FUNCT__ 
-#define __FUNCT__ "Beam::Du"
-void  Beam::Du(_p_Vec*,void*,int,int){
-	throw ErrorException(__FUNCT__," not supported yet!");
-}
-#undef __FUNCT__ 
-#define __FUNCT__ "Beam::Gradj"
-void  Beam::Gradj(_p_Vec*, void*, int, int,char*){
-	throw ErrorException(__FUNCT__," not supported yet!");
-}
-#undef __FUNCT__ 
-#define __FUNCT__ "Beam::GradjDrag"
-void  Beam::GradjDrag(_p_Vec*, void*, int,int ){
-	throw ErrorException(__FUNCT__," not supported yet!");
-}
-#undef __FUNCT__ 
-#define __FUNCT__ "Beam::GradjB"
-void  Beam::GradjB(_p_Vec*, void*, int, int){
-	throw ErrorException(__FUNCT__," not supported yet!");
-}
-#undef __FUNCT__ 
-#define __FUNCT__ "Beam::Misfit"
-double Beam::Misfit(void*,int,int ){
-	throw ErrorException(__FUNCT__," not supported yet!");
-}
-
-#undef __FUNCT__ 
-#define __FUNCT__ "Beam::GetThicknessList"
-void  Beam::GetThicknessList(double* thickness_list){
-	throw ErrorException(__FUNCT__," not supported yet!");
-}
-
-#undef __FUNCT__ 
-#define __FUNCT__ "Beam::GetParameterValue"
-void Beam::GetParameterValue(double* pvalue, double* value_list,double gauss_coord){
-
-	double l1l2[2];
-	
-	GetNodalFunctions(&l1l2[0],gauss_coord);
-
-	*pvalue=l1l2[0]*value_list[0]+l1l2[1]*value_list[1];
-}
-
-
-#undef __FUNCT__ 
-#define __FUNCT__ "Beam::GetNodalFunctions"
-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]=.5*gauss_coord+.5;
-
-	/*Second nodal function: */
-	l1l2[1]=-.5*gauss_coord+.5;
-}
-
-#undef __FUNCT__ 
-#define __FUNCT__ "Beam::GetJacobianDeterminant"
-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){
-		throw ErrorException(__FUNCT__," negative jacobian determinant!");
-	}
-	
-	*pJdet=Jdet;
-}
-
-#undef __FUNCT__ 
-#define __FUNCT__ "Beam::ComputePressure"
-void  Beam::ComputePressure(Vec p_g){
-
-	int i;
-	const int numgrids=2;
-	int doflist[numgrids];
-	double pressure[numgrids];
-	double rho_ice,g;
-	double xyz_list[numgrids][3];
-
-	/*Get node data: */
-	GetElementNodeData( &xyz_list[0][0], nodes, numgrids);
-        
-	/*Get dof list on which we will plug the pressure values: */
-	GetDofList1(&doflist[0]);
-
-	/*pressure is lithostatic: */
-	rho_ice=matpar->GetRhoIce();
-	g=matpar->GetG();
-	for(i=0;i<numgrids;i++){
-		pressure[i]=rho_ice*g*(s[i]-xyz_list[i][2]);
-	}
-	
-	/*plug local pressure values into global pressure vector: */
-	VecSetValues(p_g,numgrids,doflist,(const double*)pressure,INSERT_VALUES);
-
-}
-#undef __FUNCT__ 
-#define __FUNCT__ "Beam::MassFlux"
-double Beam::MassFlux( double* segment,double* ug){
-	throw ErrorException(__FUNCT__," not supported yet!");
-}
+/*}}}*/
Index: /issm/trunk/src/c/objects/Friction.cpp
===================================================================
--- /issm/trunk/src/c/objects/Friction.cpp	(revision 2906)
+++ /issm/trunk/src/c/objects/Friction.cpp	(revision 2907)
@@ -11,4 +11,5 @@
   --------------------------------------------------*/
 
+/*FUNCTION NewFriction{{{1*/
 Friction* NewFriction(void)
 {
@@ -17,5 +18,6 @@
 	return (Friction*)xmalloc(sizeof(Friction));
 }
-	
+/*}}}*/
+/*FUNCTION FrictionEcho{{{1*/
 void  FrictionEcho(Friction* friction){
 
@@ -49,6 +51,6 @@
 	printf("\n");
 }
-
-
+/*}}}*/
+/*FUNCTION FrictionInit {{{1*/
 /*--------------------------------------------------
 	FrictionInit
@@ -71,5 +73,6 @@
 	return 1;
 }
-
+/*}}}*/
+/*FUNCTION DeleteFriction {{{1*/
 /*--------------------------------------------------
 	DeleteFriction
@@ -86,6 +89,6 @@
 	xfree((void**)pfriction);
 }
-
-
+/*}}}*/
+/*FUNCTION FrictionGetAlpha2 {{{1*/
 /*--------------------------------------------------
 	FrictionGetAlpha2
@@ -138,5 +141,6 @@
 	}
 }
-
+/*}}}*/
+/*FUNCTION FrictionGetAlphaComplement {{{1*/
 /*--------------------------------------------------
 	FrictionGetAlphaComplement
@@ -180,2 +184,3 @@
 	}
 }
+/*}}}*/
Index: /issm/trunk/src/c/objects/Icefront.cpp
===================================================================
--- /issm/trunk/src/c/objects/Icefront.cpp	(revision 2906)
+++ /issm/trunk/src/c/objects/Icefront.cpp	(revision 2907)
@@ -19,8 +19,11 @@
 
 		
+/*Object constructors and destructor*/
+/*FUNCTION Icefront constructor {{{1*/
 Icefront::Icefront(){
 	return;
 }
-
+/*}}}*/
+/*FUNCTION Icefront creation {{{1*/
 Icefront::Icefront(char icefront_type[ICEFRONTSTRING],int icefront_sid, int icefront_mparid, int icefront_eid, int icefront_element_type, 
 		int icefront_node_ids[MAX_ICEFRONT_GRIDS],double icefront_h[MAX_ICEFRONT_GRIDS],double	icefront_b[MAX_ICEFRONT_GRIDS]){
@@ -51,85 +54,50 @@
 	return;
 }
-
+/*}}}*/
+/*FUNCTION Icefront destructor {{{1*/
 Icefront::~Icefront(){
 	return;
 }
-		
-void Icefront::Echo(void){
+/*}}}*/
+		
+/*Object marshall*/
+/*FUNCTION Icefront Demarshall {{{1*/
+void  Icefront::Demarshall(char** pmarshalled_dataset){
 
 	int i;
-	
-	printf("Icefront:\n");
-	printf("   type: %s\n",type);
-	printf("   sid: %i\n",sid);
-	
-	printf("   mparid: %i\n",mparid);
-	printf("   matpar_offset: %i\n",matpar_offset);
-	printf("   matpar: \n");
-	
-	printf("   element_type: %i\n",element_type);
-	printf("   eid: %i\n",eid);
-	printf("   element_offset: %i\n",eid);
-	printf("   element\n");
-		
-	if (strcmp(type,"segment")==0){
-		printf("   node_ids=[%i,%i]\n",node_ids[0],node_ids[1]);
-		printf("   node_offsets=[%i,%i]\n",node_offsets[0],node_offsets[1]);
-		printf("   h=[%g,%g]\n",h[0],h[1]);
-		printf("   b=[%g,%g]\n",b[0],b[1]);
-		for(i=0;i<2;i++){
-			if(nodes[i])nodes[i]->Echo();
-		}
-	}
-	else{
-		printf("   node_ids=[%i,%i,%i,%i]\n",node_ids[0],node_ids[1],node_ids[2],node_ids[3]);
-		printf("   node_offsets=[%i,%i]\n",node_offsets[0],node_offsets[1],node_offsets[2],node_offsets[3]);
-		printf("   h=[%g,%g,%g,%g]\n",h[0],h[1],h[2],h[3]);
-		printf("   b=[%g,%g,%g,%g]\n",b[0],b[1],b[2],b[3]);
-		for(i=0;i<4;i++){
-			if(nodes[i])nodes[i]->Echo();
-		}
-	}
-
+	char* marshalled_dataset=NULL;
+
+	/*recover marshalled_dataset: */
+	marshalled_dataset=*pmarshalled_dataset;
+
+	/*this time, no need to get enum type, the pointer directly points to the beginning of the 
+	 *object data (thanks to DataSet::Demarshall):*/
+
+	memcpy(&type,marshalled_dataset,sizeof(type));marshalled_dataset+=sizeof(type);
+	memcpy(&sid,marshalled_dataset,sizeof(sid));marshalled_dataset+=sizeof(sid);
+	
+	memcpy(&mparid,marshalled_dataset,sizeof(mparid));marshalled_dataset+=sizeof(mparid);
+	memcpy(&matpar_offset,marshalled_dataset,sizeof(matpar_offset));marshalled_dataset+=sizeof(matpar_offset);
+
+	memcpy(&element_type,marshalled_dataset,sizeof(element_type));marshalled_dataset+=sizeof(element_type);
+	memcpy(&eid,marshalled_dataset,sizeof(eid));marshalled_dataset+=sizeof(eid);
+	memcpy(&element_offset,marshalled_dataset,sizeof(element_offset));marshalled_dataset+=sizeof(element_offset);
+
+	memcpy(&node_ids,marshalled_dataset,sizeof(node_ids));marshalled_dataset+=sizeof(node_ids);
+	memcpy(&node_offsets,marshalled_dataset,sizeof(node_offsets));marshalled_dataset+=sizeof(node_offsets);
+	
+	memcpy(&h,marshalled_dataset,sizeof(h));marshalled_dataset+=sizeof(h);
+	memcpy(&b,marshalled_dataset,sizeof(b));marshalled_dataset+=sizeof(b);
+
+	for(i=0;i<MAX_ICEFRONT_GRIDS;i++)nodes[i]=NULL;
+	matpar=NULL;
+	element=NULL;
+
+	/*return: */
+	*pmarshalled_dataset=marshalled_dataset;
 	return;
 }
-void Icefront::DeepEcho(void){
-
-	int i;
-	
-	printf("Icefront:\n");
-	printf("   type: %s\n",type);
-	printf("   sid: %i\n",sid);
-	
-	printf("   mparid: %i\n",mparid);
-	printf("   matpar_offset: %i\n",matpar_offset);
-	printf("   matpar: \n");
-	
-	printf("   element_type: %i\n",element_type);
-	printf("   eid: %i\n",eid);
-	printf("   element_offset: %i\n",eid);
-	printf("   element\n");
-		
-	if (strcmp(type,"segment")==0){
-		printf("   node_ids=[%i,%i]\n",node_ids[0],node_ids[1]);
-		printf("   node_offsets=[%i,%i]\n",node_offsets[0],node_offsets[1]);
-		printf("   h=[%g,%g]\n",h[0],h[1]);
-		printf("   b=[%g,%g]\n",b[0],b[1]);
-		for(i=0;i<2;i++){
-			if(nodes[i])nodes[i]->Echo();
-		}
-	}
-	else{
-		printf("   node_ids=[%i,%i,%i,%i]\n",node_ids[0],node_ids[1],node_ids[2],node_ids[3]);
-		printf("   node_offsets=[%i,%i]\n",node_offsets[0],node_offsets[1],node_offsets[2],node_offsets[3]);
-		printf("   h=[%g,%g,%g,%g]\n",h[0],h[1],h[2],h[3]);
-		printf("   b=[%g,%g,%g,%g]\n",b[0],b[1],b[2],b[3]);
-		for(i=0;i<4;i++){
-			if(nodes[i])nodes[i]->Echo();
-		}
-	}
-
-	return;
-}		
+/*}}}*/
+/*FUNCTION Icefront Marshall {{{1*/
 void  Icefront::Marshall(char** pmarshalled_dataset){
 
@@ -166,5 +134,6 @@
 	return;
 }
-		
+/*}}}*/
+/*FUNCTION Icefront MarshallSize{{{1*/
 int   Icefront::MarshallSize(){
 
@@ -182,60 +151,8 @@
 		sizeof(int); //sizeof(int) for enum type
 }
-
-char* Icefront::GetName(void){
-	return "icefront";
-}
-		
-
-void  Icefront::Demarshall(char** pmarshalled_dataset){
-
-	int i;
-	char* marshalled_dataset=NULL;
-
-	/*recover marshalled_dataset: */
-	marshalled_dataset=*pmarshalled_dataset;
-
-	/*this time, no need to get enum type, the pointer directly points to the beginning of the 
-	 *object data (thanks to DataSet::Demarshall):*/
-
-	memcpy(&type,marshalled_dataset,sizeof(type));marshalled_dataset+=sizeof(type);
-	memcpy(&sid,marshalled_dataset,sizeof(sid));marshalled_dataset+=sizeof(sid);
-	
-	memcpy(&mparid,marshalled_dataset,sizeof(mparid));marshalled_dataset+=sizeof(mparid);
-	memcpy(&matpar_offset,marshalled_dataset,sizeof(matpar_offset));marshalled_dataset+=sizeof(matpar_offset);
-
-	memcpy(&element_type,marshalled_dataset,sizeof(element_type));marshalled_dataset+=sizeof(element_type);
-	memcpy(&eid,marshalled_dataset,sizeof(eid));marshalled_dataset+=sizeof(eid);
-	memcpy(&element_offset,marshalled_dataset,sizeof(element_offset));marshalled_dataset+=sizeof(element_offset);
-
-	memcpy(&node_ids,marshalled_dataset,sizeof(node_ids));marshalled_dataset+=sizeof(node_ids);
-	memcpy(&node_offsets,marshalled_dataset,sizeof(node_offsets));marshalled_dataset+=sizeof(node_offsets);
-	
-	memcpy(&h,marshalled_dataset,sizeof(h));marshalled_dataset+=sizeof(h);
-	memcpy(&b,marshalled_dataset,sizeof(b));marshalled_dataset+=sizeof(b);
-
-	for(i=0;i<MAX_ICEFRONT_GRIDS;i++)nodes[i]=NULL;
-	matpar=NULL;
-	element=NULL;
-
-	/*return: */
-	*pmarshalled_dataset=marshalled_dataset;
-	return;
-}
-
-int Icefront::Enum(void){
-
-	return IcefrontEnum();
-
-}
-
-int    Icefront::GetId(void){ return sid; }
-
-int    Icefront::MyRank(void){ 
-	extern int my_rank;
-	return my_rank; 
-}
-void  Icefront::DistributeNumDofs(int* numdofspernode,int analysis_type,int sub_analysis_type){return;}
-		
+/*}}}*/
+
+/*Object functions*/
+/*FUNCTION Icefront Configure {{{1*/
 #undef __FUNCT__ 
 #define __FUNCT__ "Icefront::Configure"
@@ -271,7 +188,11 @@
 
 }
-
-
-
+/*}}}*/
+/*FUNCTION Icefront copy {{{1*/
+Object* Icefront::copy() {
+	return new Icefront(*this); 
+}
+/*}}}*/
+/*FUNCTION Icefront CreateKMatrix {{{1*/
 #undef __FUNCT__ 
 #define __FUNCT__ "Icefront::CreateKMatrix"
@@ -284,5 +205,6 @@
 
 }
-
+/*}}}*/
+/*FUNCTION Icefront CreatePVector {{{1*/
 #undef __FUNCT__ 
 #define __FUNCT__ "Icefront::CreatePVector"
@@ -314,6 +236,6 @@
 	}
 }
-
-
+/*}}}*/
+/*FUNCTION Icefront CreatePVectorDiagnosticHoriz {{{1*/
 #undef __FUNCT__ 
 #define __FUNCT__ "Icefront::CreatePVectorDiagnosticHoriz"
@@ -328,5 +250,6 @@
 	}
 }	
-
+/*}}}*/
+/*FUNCTION Icefront CreatePVectorDiagnosticHorizStokes {{{1*/
 #undef __FUNCT__ 
 #define __FUNCT__ "Icefront::CreatePVectorDiagnosticHorizSegment"
@@ -443,5 +366,6 @@
 
 }
-
+/*}}}*/
+/*FUNCTION Icefront CreatePVectorDiagnosticHorizQuad {{{1*/
 #undef __FUNCT__ 
 #define __FUNCT__ "Icefont::CreatePVectorDiagnosticHorizQuad"
@@ -605,6 +529,6 @@
 
 }
-
-
+/*}}}*/
+/*FUNCTION Icefront CreatePVectorDiagnosticStokes {{{1*/
 #undef __FUNCT__ 
 #define __FUNCT__ "Icefont::CreatePVectorDiagnosticStokes"
@@ -767,34 +691,97 @@
 
 }
-
-
-#undef __FUNCT__ 
-#define __FUNCT__ "Icefront::UpdateFromInputs"
-void  Icefront::UpdateFromInputs(void* vinputs){
+/*}}}*/
+/*FUNCTION Icefront DeepEcho {{{1*/
+void Icefront::DeepEcho(void){
 
 	int i;
-	int numberofdofspernode;
-
-	/*element: */
-	double* h_param=NULL;
-	double* b_param=NULL;
-	int  dofs[1]={0};
-
-	ParameterInputs* inputs=NULL;	
-
-	inputs=(ParameterInputs*)vinputs;
-
-	if(strcmp(type,"quad")==0){
-		inputs->Recover("thickness",&h[0],1,dofs,4,(void**)nodes);
-		inputs->Recover("bed",&b[0],1,dofs,4,(void**)nodes);
+	
+	printf("Icefront:\n");
+	printf("   type: %s\n",type);
+	printf("   sid: %i\n",sid);
+	
+	printf("   mparid: %i\n",mparid);
+	printf("   matpar_offset: %i\n",matpar_offset);
+	printf("   matpar: \n");
+	
+	printf("   element_type: %i\n",element_type);
+	printf("   eid: %i\n",eid);
+	printf("   element_offset: %i\n",eid);
+	printf("   element\n");
+		
+	if (strcmp(type,"segment")==0){
+		printf("   node_ids=[%i,%i]\n",node_ids[0],node_ids[1]);
+		printf("   node_offsets=[%i,%i]\n",node_offsets[0],node_offsets[1]);
+		printf("   h=[%g,%g]\n",h[0],h[1]);
+		printf("   b=[%g,%g]\n",b[0],b[1]);
+		for(i=0;i<2;i++){
+			if(nodes[i])nodes[i]->Echo();
+		}
 	}
 	else{
-		inputs->Recover("thickness",&h[0],1,dofs,2,(void**)nodes);
-		inputs->Recover("bed",&h[0],1,dofs,2,(void**)nodes);
-	}
-
-
-}
-
+		printf("   node_ids=[%i,%i,%i,%i]\n",node_ids[0],node_ids[1],node_ids[2],node_ids[3]);
+		printf("   node_offsets=[%i,%i]\n",node_offsets[0],node_offsets[1],node_offsets[2],node_offsets[3]);
+		printf("   h=[%g,%g,%g,%g]\n",h[0],h[1],h[2],h[3]);
+		printf("   b=[%g,%g,%g,%g]\n",b[0],b[1],b[2],b[3]);
+		for(i=0;i<4;i++){
+			if(nodes[i])nodes[i]->Echo();
+		}
+	}
+
+	return;
+}		
+/*}}}*/
+/*FUNCTION Icefront DistributeNumDofs {{{1*/
+void  Icefront::DistributeNumDofs(int* numdofspernode,int analysis_type,int sub_analysis_type){return;}
+		
+/*}}}*/
+/*FUNCTION Icefront Echo {{{1*/
+void Icefront::Echo(void){
+
+	int i;
+	
+	printf("Icefront:\n");
+	printf("   type: %s\n",type);
+	printf("   sid: %i\n",sid);
+	
+	printf("   mparid: %i\n",mparid);
+	printf("   matpar_offset: %i\n",matpar_offset);
+	printf("   matpar: \n");
+	
+	printf("   element_type: %i\n",element_type);
+	printf("   eid: %i\n",eid);
+	printf("   element_offset: %i\n",eid);
+	printf("   element\n");
+		
+	if (strcmp(type,"segment")==0){
+		printf("   node_ids=[%i,%i]\n",node_ids[0],node_ids[1]);
+		printf("   node_offsets=[%i,%i]\n",node_offsets[0],node_offsets[1]);
+		printf("   h=[%g,%g]\n",h[0],h[1]);
+		printf("   b=[%g,%g]\n",b[0],b[1]);
+		for(i=0;i<2;i++){
+			if(nodes[i])nodes[i]->Echo();
+		}
+	}
+	else{
+		printf("   node_ids=[%i,%i,%i,%i]\n",node_ids[0],node_ids[1],node_ids[2],node_ids[3]);
+		printf("   node_offsets=[%i,%i]\n",node_offsets[0],node_offsets[1],node_offsets[2],node_offsets[3]);
+		printf("   h=[%g,%g,%g,%g]\n",h[0],h[1],h[2],h[3]);
+		printf("   b=[%g,%g,%g,%g]\n",b[0],b[1],b[2],b[3]);
+		for(i=0;i<4;i++){
+			if(nodes[i])nodes[i]->Echo();
+		}
+	}
+
+	return;
+}
+/*}}}*/
+/*FUNCTION Icefront Enum {{{1*/
+int Icefront::Enum(void){
+
+	return IcefrontEnum();
+
+}
+/*}}}*/
+/*FUNCTION Icefront GetDofList{{{1*/
 #undef __FUNCT__ 
 #define __FUNCT__ "Icefront::GetDofList"
@@ -828,76 +815,20 @@
 
 }
-
-
-#undef __FUNCT__ 
-#define __FUNCT__ "Icefront::SegmentPressureLoad"
-
-void Icefront::SegmentPressureLoad(double* pe_g,double rho_water,double rho_ice,double gravity, double* thickness_list, double* bed_list, double* normal,double length,int fill){
-
-	double   nx,ny;
-	double   h1,h2,b1,b2;
-	int      num_gauss;
-	double*  segment_gauss_coord=NULL;
-	double*  gauss_weights=NULL;
-	int      ig;
-	double   Jdet;
-	double   thickness,bed;
-	double   ice_pressure,water_pressure,air_pressure;
-	double   surface_under_water,base_under_water;
-	double   pressure;
-
-	nx=normal[0];
-	ny=normal[1];
-
-	//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);
-
-	//recover thickness and bed at two grids
-	h1=thickness_list[0];
-	h2=thickness_list[1];
-	b1=bed_list[0];
-	b2=bed_list[1];
-
-	//compute Jacobian of segment
-	Jdet=1./2.*length;
-
-	for(ig=0;ig<num_gauss;ig++){
-
-		thickness=h1*(1+segment_gauss_coord[ig])/2+h2*(1-segment_gauss_coord[ig])/2;
-		bed=b1*(1+segment_gauss_coord[ig])/2+b2*(1-segment_gauss_coord[ig])/2;
-
-		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{
-			throw ErrorException(__FUNCT__," fill type not supported yet");
-		}
-
-		pressure = ice_pressure + water_pressure + air_pressure;
-
-		pe_g[2*0+0]+= pressure*Jdet*gauss_weights[ig]*(1+segment_gauss_coord[ig])/2*nx;
-		pe_g[2*0+1]+= pressure*Jdet*gauss_weights[ig]*(1+segment_gauss_coord[ig])/2*ny;
-		pe_g[2*1+0]+= pressure*Jdet*gauss_weights[ig]*(1-segment_gauss_coord[ig])/2*nx;
-		pe_g[2*1+1]+= pressure*Jdet*gauss_weights[ig]*(1-segment_gauss_coord[ig])/2*ny;
-
-	} //for(ig=0;ig<num_gauss;ig++)
-
-	xfree((void**)&segment_gauss_coord);
-	xfree((void**)&gauss_weights);
-}
-
+/*}}}*/
+/*FUNCTION Icefront GetId {{{1*/
+int    Icefront::GetId(void){ return sid; }
+/*}}}*/
+/*FUNCTION Icefront GetName {{{1*/
+char* Icefront::GetName(void){
+	return "icefront";
+}
+/*}}}*/
+/*FUNCTION Icefront MyRank {{{1*/
+int    Icefront::MyRank(void){ 
+	extern int my_rank;
+	return my_rank; 
+}
+/*}}}*/
+/*FUNCTION Icefront PenaltyCreateKMatrix {{{1*/
 #undef __FUNCT__ 
 #define __FUNCT__ "Icefront::PenaltyCreateKMatrix"
@@ -905,5 +836,6 @@
 	/*do nothing: */
 }
-		
+/*}}}*/
+/*FUNCTION Icefront PenaltyCreatePVector{{{1*/
 #undef __FUNCT__ 
 #define __FUNCT__ "Icefront::PenaltyCreatePVector"
@@ -911,10 +843,6 @@
 	/*do nothing: */
 }
-
-Object* Icefront::copy() {
-	return new Icefront(*this); 
-}
-
-
+/*}}}*/
+/*FUNCTION Icefront QuadPressureLoad {{{1*/
 #undef __FUNCT__ 
 #define __FUNCT__ "Icefront::QuadPressureLoad"
@@ -1173,5 +1101,6 @@
 	delete tria;
 }
-
+/*}}}*/
+/*FUNCTION Icefront QuadPressureLoadStokes {{{1*/
 #undef __FUNCT__ 
 #define __FUNCT__ "Icefront::QuadPressureLoadStokes"
@@ -1430,2 +1359,103 @@
 	delete tria;
 }
+/*}}}*/
+/*FUNCTION Icefront SegmentPressureLoad {{{1*/
+#undef __FUNCT__ 
+#define __FUNCT__ "Icefront::SegmentPressureLoad"
+
+void Icefront::SegmentPressureLoad(double* pe_g,double rho_water,double rho_ice,double gravity, double* thickness_list, double* bed_list, double* normal,double length,int fill){
+
+	double   nx,ny;
+	double   h1,h2,b1,b2;
+	int      num_gauss;
+	double*  segment_gauss_coord=NULL;
+	double*  gauss_weights=NULL;
+	int      ig;
+	double   Jdet;
+	double   thickness,bed;
+	double   ice_pressure,water_pressure,air_pressure;
+	double   surface_under_water,base_under_water;
+	double   pressure;
+
+	nx=normal[0];
+	ny=normal[1];
+
+	//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);
+
+	//recover thickness and bed at two grids
+	h1=thickness_list[0];
+	h2=thickness_list[1];
+	b1=bed_list[0];
+	b2=bed_list[1];
+
+	//compute Jacobian of segment
+	Jdet=1./2.*length;
+
+	for(ig=0;ig<num_gauss;ig++){
+
+		thickness=h1*(1+segment_gauss_coord[ig])/2+h2*(1-segment_gauss_coord[ig])/2;
+		bed=b1*(1+segment_gauss_coord[ig])/2+b2*(1-segment_gauss_coord[ig])/2;
+
+		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{
+			throw ErrorException(__FUNCT__," fill type not supported yet");
+		}
+
+		pressure = ice_pressure + water_pressure + air_pressure;
+
+		pe_g[2*0+0]+= pressure*Jdet*gauss_weights[ig]*(1+segment_gauss_coord[ig])/2*nx;
+		pe_g[2*0+1]+= pressure*Jdet*gauss_weights[ig]*(1+segment_gauss_coord[ig])/2*ny;
+		pe_g[2*1+0]+= pressure*Jdet*gauss_weights[ig]*(1-segment_gauss_coord[ig])/2*nx;
+		pe_g[2*1+1]+= pressure*Jdet*gauss_weights[ig]*(1-segment_gauss_coord[ig])/2*ny;
+
+	} //for(ig=0;ig<num_gauss;ig++)
+
+	xfree((void**)&segment_gauss_coord);
+	xfree((void**)&gauss_weights);
+}
+/*}}}*/
+/*FUNCTION Icefront UpdateFromInputs {{{1*/
+#undef __FUNCT__ 
+#define __FUNCT__ "Icefront::UpdateFromInputs"
+void  Icefront::UpdateFromInputs(void* vinputs){
+
+	int i;
+	int numberofdofspernode;
+
+	/*element: */
+	double* h_param=NULL;
+	double* b_param=NULL;
+	int  dofs[1]={0};
+
+	ParameterInputs* inputs=NULL;	
+
+	inputs=(ParameterInputs*)vinputs;
+
+	if(strcmp(type,"quad")==0){
+		inputs->Recover("thickness",&h[0],1,dofs,4,(void**)nodes);
+		inputs->Recover("bed",&b[0],1,dofs,4,(void**)nodes);
+	}
+	else{
+		inputs->Recover("thickness",&h[0],1,dofs,2,(void**)nodes);
+		inputs->Recover("bed",&h[0],1,dofs,2,(void**)nodes);
+	}
+
+
+}
+/*}}}*/
Index: /issm/trunk/src/c/objects/Input.cpp
===================================================================
--- /issm/trunk/src/c/objects/Input.cpp	(revision 2906)
+++ /issm/trunk/src/c/objects/Input.cpp	(revision 2907)
@@ -19,8 +19,11 @@
 
 		
+/*Object constructors and destructor*/
+/*FUNCTION Input::constructor {{{1*/
 Input::Input(){
 	return;
 }
-
+/*}}}*/
+/*FUNCTION Input::destructor {{{1*/
 Input::~Input(){
 
@@ -28,5 +31,32 @@
 	return;
 }
-
+/*}}}*/
+/*FUNCTION Input::creation(char* input_name,int input_integer) {{{1*/
+Input::Input(char* input_name,int input_integer){
+
+	strcpy(name,input_name);
+	type=INTEGER;
+	integer=input_integer;
+
+	ndof=0;
+	numberofnodes=0;
+	vector=NULL;
+
+ar* input_name,double input_scalar}
+/*}}}*/
+/*FUNCTION Input::creation(char* input_name,double input_scalar) {{{1*/
+Input::Input(char* input_name,double input_scalar){
+
+	strcpy(name,input_name);
+	type=DOUBLE;
+	scalar=input_scalar;
+
+	ndof=0;
+	numberofnodes=0;
+	vector=NULL;
+
+}
+/*}}}*/
+/*FUNCTION Input::creation(char* input_name,char* input_string) {{{1*/
 Input::Input(char* input_name,char* input_string){
 	
@@ -38,29 +68,6 @@
 	vector=NULL;
 }
-		
-Input::Input(char* input_name,int input_integer){
-
-	strcpy(name,input_name);
-	type=INTEGER;
-	integer=input_integer;
-
-	ndof=0;
-	numberofnodes=0;
-	vector=NULL;
-
-}
-		
-Input::Input(char* input_name,double input_scalar){
-
-	strcpy(name,input_name);
-	type=DOUBLE;
-	scalar=input_scalar;
-
-	ndof=0;
-	numberofnodes=0;
-	vector=NULL;
-
-}
-
+/*}}}*/
+/*FUNCTION Input::creation(char* input_name,double* input_vector,int input_ndof,int input_numberofnodes) {{{1*/
 Input::Input(char* input_name,double* input_vector,int input_ndof,int input_numberofnodes){
 
@@ -73,5 +80,6 @@
 
 }
-
+/*}}}*/
+/*FUNCTION Input::creation(char* input_name,Vec input_vector,int input_ndof, int input_numberofnode) {{{1*/
 Input::Input(char* input_name,Vec input_vector,int input_ndof, int input_numberofnodes){
 
@@ -83,18 +91,38 @@
 
 }
-
-void Input::Echo(void){
-
-	int i,j;
-	
-	printf("Input:\n");
-	printf("   name: %s\n",name);
-	printf("   type: %i\n",type);
-	if(type==STRING) printf("   string: %s\n",string);
-	if(type==INTEGER) printf("   integer: %i\n",integer);
-	if(type==DOUBLE) printf("   double: %g\n",scalar);
-	if(type==DOUBLEVEC) printf("   doublevec:%p\n",vector);
-}
-
+/*}}}*/
+
+/*Object marshall*/
+/*FUNCTION Input::Demarshall {{{1*/
+#undef __FUNCT__
+#define __FUNCT__ "Input::Demarshall"
+void  Input::Demarshall(char** pmarshalled_dataset){
+	throw ErrorException(__FUNCT__,"not supported yet!");
+}
+/*}}}*/
+/*FUNCTION Input::Marshall {{{1*/
+#undef __FUNCT__
+#define __FUNCT__ "Input::Marshall"
+void  Input::Marshall(char** pmarshalled_dataset){
+	throw ErrorException(__FUNCT__,"not supported yet!");
+}
+/*}}}*/
+/*FUNCTION Input::MarshallSize {{{1*/
+#undef __FUNCT__
+#define __FUNCT__ "Input::MarshallSize"
+int   Input::MarshallSize(){
+	throw ErrorException(__FUNCT__,"not supported yet!");
+}
+/*}}}*/
+
+/*Object functions*/
+/*FUNCTION Input::copy {{{1*/
+Object* Input::copy() {
+	
+	return new Input(*this); 
+
+}
+/*}}}*/
+/*FUNCTION Input::DeepEcho {{{1*/
 void Input::DeepEcho(void){
 
@@ -114,4 +142,20 @@
 	}
 }		
+/*}}}*/
+/*FUNCTION Input::Echo {{{1*/
+void Input::Echo(void){
+
+	int i,j;
+	
+	printf("Input:\n");
+	printf("   name: %s\n",name);
+	printf("   type: %i\n",type);
+	if(type==STRING) printf("   string: %s\n",string);
+	if(type==INTEGER) printf("   integer: %i\n",integer);
+	if(type==DOUBLE) printf("   double: %g\n",scalar);
+	if(type==DOUBLEVEC) printf("   doublevec:%p\n",vector);
+}
+/*}}}*/
+/*FUNCTION Input::Enum {{{1*/
 int Input::Enum(void){
 
@@ -119,118 +163,6 @@
 
 }
-
-int    Input::MyRank(void){ 
-	extern int my_rank;
-	return my_rank; 
-}
-		
-
-int   Input::IsSameName(char* input_name){
-	if (strcmp(input_name,name)==0)return 1;
-	else return 0;
-}
-		
-
-void  Input::Recover(double* valuesin, int ndofin, int* dofsin,int numnodesin,void** vpnodesin){
-
-	//ex: 
-	// double values[6];
-	// int ndof=2;
-	// int dofs[2]={0,2}    //dofs 0 and 2
-	// for a tria, 3 nodes
-	//input->Recover(&values[0],ndof,dofs,3,this->nodes); //ie: recover values for dofs 0 and 2 of all nodes of the tria
-
-	int* doflist=NULL;
-	int  dof1;
-	int  dof;
-	int  i,j;
-	
-	Node** nodes=NULL;
-	Node* node=NULL;
-
-	/*recover pointer to nodes: */
-	nodes=(Node**)vpnodesin;
-	
-	if (type!=DOUBLEVEC) throw ErrorException(__FUNCT__,exprintf("%s%i%s"," cannot recover values from a ",type," input type"));
-	
-	/*Ok, we are trying to fill values. values is of size ndof*numnodes.  The dofs 
-	 * for picking up the values are found in the numnodes nodes. We firt build the dof list: */
-	doflist=(int*)xmalloc(ndofin*numnodesin*sizeof(int));
-
-	for(i=0;i<numnodesin;i++){
-		node=nodes[i];
-
-		/*Get 1d dof of this node: */
-		dof1=node->GetDofList1();
-
-		for(j=0;j<ndofin;j++){
-			dof=dofsin[j];
-			doflist[i*ndofin+j]=dof1*ndof+dof;
-		}
-	}
-	
-	for(i=0;i<ndofin*numnodesin;i++){
-		valuesin[i]=vector[doflist[i]];
-	}
-
-	xfree((void**)&doflist);
-}
-
-Object* Input::copy() {
-	
-	return new Input(*this); 
-
-}
-
-
-void  Input::Recover(int* pinteger){
-	*pinteger=integer;
-}
-	
-void  Input::Recover(char** pstring){
-	char* outstring=NULL;
-
-	outstring=(char*)xmalloc((strlen(string)+1)*sizeof(char));
-	strcpy(outstring,string);
-	*pstring=outstring;
-
-}
-
-void  Input::Recover(double* pdouble){
-
-	*pdouble=scalar;
-
-}
-
-#undef __FUNCT__
-#define __FUNCT__ "Input::GetId"
-int   Input::GetId(){
-	throw ErrorException(__FUNCT__,"not supported yet!");
-}
-
-#undef __FUNCT__
-#define __FUNCT__ "Input::Marshall"
-void  Input::Marshall(char** pmarshalled_dataset){
-	throw ErrorException(__FUNCT__,"not supported yet!");
-}
-
-#undef __FUNCT__
-#define __FUNCT__ "Input::MarshallSize"
-int   Input::MarshallSize(){
-	throw ErrorException(__FUNCT__,"not supported yet!");
-}
-
-#undef __FUNCT__
-#define __FUNCT__ "Input::GetName"
-char* Input::GetName(){
-	throw ErrorException(__FUNCT__,"not supported yet!");
-}
-		
-#undef __FUNCT__
-#define __FUNCT__ "Input::Demarshall"
-void  Input::Demarshall(char** pmarshalled_dataset){
-	throw ErrorException(__FUNCT__,"not supported yet!");
-}
-
+/*}}}*/
+/*FUNCTION Input::Get {{{1*/
 #undef __FUNCT__
 #define __FUNCT__ "Input::Get"
@@ -279,2 +211,102 @@
 	return outvector;
 }
+/*}}}*/
+/*FUNCTION Input::GetId {{{1*/
+#undef __FUNCT__
+#define __FUNCT__ "Input::GetId"
+int   Input::GetId(){
+	throw ErrorException(__FUNCT__,"not supported yet!");
+}
+/*}}}*/
+/*FUNCTION Input::GetName {{{1*/
+#undef __FUNCT__
+#define __FUNCT__ "Input::GetName"
+char* Input::GetName(){
+	throw ErrorException(__FUNCT__,"not supported yet!");
+}
+/*}}}*/
+/*FUNCTION Input::IsSameName {{{1*/
+int   Input::IsSameName(char* input_name){
+	if (strcmp(input_name,name)==0)return 1;
+	else return 0;
+}
+/*}}}*/
+/*FUNCTION Input::MyRank {{{1*/
+int    Input::MyRank(void){ 
+	extern int my_rank;
+	return my_rank; 
+}
+/*}}}*/
+/*FUNCTION Input::Recover(double* valuesin, int ndofin, int* dofsin,int numnodesin,void** vpnodesin) {{{1*/
+void  Input::Recover(double* valuesin, int ndofin, int* dofsin,int numnodesin,void** vpnodesin){
+
+	//ex: 
+	// double values[6];
+	// int ndof=2;
+	// int dofs[2]={0,2}    //dofs 0 and 2
+	// for a tria, 3 nodes
+	//input->Recover(&values[0],ndof,dofs,3,this->nodes); //ie: recover values for dofs 0 and 2 of all nodes of the tria
+
+	int* doflist=NULL;
+	int  dof1;
+	int  dof;
+	int  i,j;
+	
+	Node** nodes=NULL;
+	Node* node=NULL;
+
+	/*recover pointer to nodes: */
+	nodes=(Node**)vpnodesin;
+	
+	if (type!=DOUBLEVEC) throw ErrorException(__FUNCT__,exprintf("%s%i%s"," cannot recover values from a ",type," input type"));
+	
+	/*Ok, we are trying to fill values. values is of size ndof*numnodes.  The dofs 
+	 * for picking up the values are found in the numnodes nodes. We firt build the dof list: */
+	doflist=(int*)xmalloc(ndofin*numnodesin*sizeof(int));
+
+	for(i=0;i<numnodesin;i++){
+		node=nodes[i];
+
+		/*Get 1d dof of this node: */
+		dof1=node->GetDofList1();
+
+		for(j=0;j<ndofin;j++){
+			dof=dofsin[j];
+			doflist[i*ndofin+j]=dof1*ndof+dof;
+		}
+	}
+	
+	for(i=0;i<ndofin*numnodesin;i++){
+		valuesin[i]=vector[doflist[i]];
+	}
+
+	xfree((void**)&doflist);
+}
+/*}}}*/
+/*FUNCTION Input::Recover(int* pinteger) {{{1*/
+void  Input::Recover(int* pinteger){
+	*pinteger=integer;
+}
+/*}}}*/
+/*FUNCTION Input::Recover(double* pdouble) {{{1*/
+void  Input::Recover(double* pdouble){
+
+	*pdouble=scalar;
+
+}
+/*}}}*/
+/*FUNCTION Input::Recover(char** pstring) {{{1*/
+void  Input::Recover(char** pstring){
+	char* outstring=NULL;
+
+	outstring=(char*)xmalloc((strlen(string)+1)*sizeof(char));
+	strcpy(outstring,string);
+	*pstring=outstring;
+
+}
+/*}}}*/
+
+
+
+		
+
Index: /issm/trunk/src/c/objects/Matice.cpp
===================================================================
--- /issm/trunk/src/c/objects/Matice.cpp	(revision 2906)
+++ /issm/trunk/src/c/objects/Matice.cpp	(revision 2907)
@@ -17,7 +17,11 @@
 
 		
+/*Object constructors and destructor*/
+/*FUNCTION Matice::constructor {{{1*/
 Matice::Matice(){
 	return;
 }
+/*}}}*/
+/*FUNCTION Matice::creation {{{1*/
 Matice::Matice(int matice_mid,double matice_B,double matice_n){
 	mid=matice_mid;
@@ -26,9 +30,71 @@
 	return;
 }
-
+/*}}}*/
+/*FUNCTION Matice::destructor {{{1*/
 Matice::~Matice(){
 	return;
 }
-void Matice::Echo(void){
+/*}}}*/
+
+/*Object marshall*/
+/*FUNCTION Matice::Marshall {{{1*/
+void  Matice::Marshall(char** pmarshalled_dataset){
+
+	char* marshalled_dataset=NULL;
+	int   enum_type=0;
+
+	/*recover marshalled_dataset: */
+	marshalled_dataset=*pmarshalled_dataset;
+
+	/*get enum type of Matice: */
+	enum_type=MaticeEnum();
+	
+	/*marshall enum: */
+	memcpy(marshalled_dataset,&enum_type,sizeof(enum_type));marshalled_dataset+=sizeof(enum_type);
+	
+	/*marshall Matice data: */
+	memcpy(marshalled_dataset,&mid,sizeof(mid));marshalled_dataset+=sizeof(mid);
+	memcpy(marshalled_dataset,&B,sizeof(B));marshalled_dataset+=sizeof(B);
+	memcpy(marshalled_dataset,&n,sizeof(n));marshalled_dataset+=sizeof(n);
+
+	*pmarshalled_dataset=marshalled_dataset;
+	return;
+}
+/*}}}*/
+/*FUNCTION Matice::MarshallSize{{{1*/
+int   Matice::MarshallSize(){
+
+	return sizeof(mid)+sizeof(B)+sizeof(n)+sizeof(int); //sizeof(int) for enum type
+}
+/*}}}*/
+/*FUNCTION Matice::Demarshall {{{1*/
+void  Matice::Demarshall(char** pmarshalled_dataset){
+
+	char* marshalled_dataset=NULL;
+
+	/*recover marshalled_dataset: */
+	marshalled_dataset=*pmarshalled_dataset;
+
+	/*this time, no need to get enum type, the pointer directly points to the beginning of the 
+	 *object data (thanks to DataSet::Demarshall):*/
+
+	memcpy(&mid,marshalled_dataset,sizeof(mid));marshalled_dataset+=sizeof(mid);
+	memcpy(&B,marshalled_dataset,sizeof(B));marshalled_dataset+=sizeof(B);
+	memcpy(&n,marshalled_dataset,sizeof(n));marshalled_dataset+=sizeof(n);
+
+	/*return: */
+	*pmarshalled_dataset=marshalled_dataset;
+	return;
+}
+/*}}}*/
+
+/*Object functions*/
+/*FUNCTION Matice::copy {{{1*/
+Object* Matice::copy() {
+	return new Matice(*this); 
+}
+/*}}}*/
+/*FUNCTION Matice::DeepEcho {{{1*/
+void Matice::DeepEcho(void){
 
 	printf("Matice:\n");
@@ -37,6 +103,11 @@
 	printf("   n: %g\n",n);
 	return;
-}
-void Matice::DeepEcho(void){
+}		
+/*}}}*/
+/*FUNCTION Matice::DistributeNumDofs {{{1*/
+void  Matice::DistributeNumDofs(int* numdofspernode,int analysis_type){return;}
+/*}}}*/
+/*FUNCTION Matice::Echo {{{1*/
+void Matice::Echo(void){
 
 	printf("Matice:\n");
@@ -45,86 +116,32 @@
 	printf("   n: %g\n",n);
 	return;
-}		
-void  Matice::Marshall(char** pmarshalled_dataset){
-
-	char* marshalled_dataset=NULL;
-	int   enum_type=0;
-
-	/*recover marshalled_dataset: */
-	marshalled_dataset=*pmarshalled_dataset;
-
-	/*get enum type of Matice: */
-	enum_type=MaticeEnum();
-	
-	/*marshall enum: */
-	memcpy(marshalled_dataset,&enum_type,sizeof(enum_type));marshalled_dataset+=sizeof(enum_type);
-	
-	/*marshall Matice data: */
-	memcpy(marshalled_dataset,&mid,sizeof(mid));marshalled_dataset+=sizeof(mid);
-	memcpy(marshalled_dataset,&B,sizeof(B));marshalled_dataset+=sizeof(B);
-	memcpy(marshalled_dataset,&n,sizeof(n));marshalled_dataset+=sizeof(n);
-
-	*pmarshalled_dataset=marshalled_dataset;
-	return;
-}
-		
-int   Matice::MarshallSize(){
-
-	return sizeof(mid)+sizeof(B)+sizeof(n)+sizeof(int); //sizeof(int) for enum type
-}
-
+}
+/*}}}*/
+/*FUNCTION Matice::Enum {{{1*/
+int Matice::Enum(void){
+
+	return MaticeEnum();
+
+}
+/*}}}*/
+/*FUNCTION Matice::GetB {{{1*/
+double Matice::GetB(){
+	return B;
+}
+/*}}}*/
+/*FUNCTION Matice::GetId {{{1*/
+int    Matice::GetId(void){ return mid; }
+/*}}}*/
+/*FUNCTION Matice::GetN {{{1*/
+double Matice::GetN(){
+	return n;
+}
+/*}}}*/
+/*FUNCTION Matice::GetName {{{1*/
 char* Matice::GetName(void){
 	return "matice";
 }
-		
-
-void  Matice::Demarshall(char** pmarshalled_dataset){
-
-	char* marshalled_dataset=NULL;
-
-	/*recover marshalled_dataset: */
-	marshalled_dataset=*pmarshalled_dataset;
-
-	/*this time, no need to get enum type, the pointer directly points to the beginning of the 
-	 *object data (thanks to DataSet::Demarshall):*/
-
-	memcpy(&mid,marshalled_dataset,sizeof(mid));marshalled_dataset+=sizeof(mid);
-	memcpy(&B,marshalled_dataset,sizeof(B));marshalled_dataset+=sizeof(B);
-	memcpy(&n,marshalled_dataset,sizeof(n));marshalled_dataset+=sizeof(n);
-
-	/*return: */
-	*pmarshalled_dataset=marshalled_dataset;
-	return;
-}
-
-int Matice::Enum(void){
-
-	return MaticeEnum();
-
-}
-
-int    Matice::GetId(void){ return mid; }
-
-int    Matice::MyRank(void){ 
-	extern int my_rank;
-	return my_rank; 
-}
-void  Matice::DistributeNumDofs(int* numdofspernode,int analysis_type){return;}
-
-#undef __FUNCT__ 
-#define __FUNCT__ "Matice::UpdateFromInputs"
-void  Matice::UpdateFromInputs(void* inputs){
-	
-	//throw ErrorException(__FUNCT__," not supported yet!");
-
-}
-		
-#undef __FUNCT__ 
-#define __FUNCT__ "Matice::SetB"
-void  Matice::SetB(double B_param){
-	B=B_param;
-}
-		
-
+/*}}}*/
+/*FUNCTION Matice::GetViscosity2d {{{1*/
 #undef __FUNCT__ 
 #define __FUNCT__ "MatIce::GetViscosity2d"
@@ -188,65 +205,6 @@
 	*pviscosity=viscosity;
 }
-
-Object* Matice::copy() {
-	return new Matice(*this); 
-}
-
-
-#undef __FUNCT__ 
-#define __FUNCT__ "MatIce::GetViscosityComplement"
-void  Matice::GetViscosityComplement(double* pviscosity_complement, double* epsilon){
-
-	/*Return viscosity accounting for steady state power law creep [Thomas and MacAyeal, 1982]: 
-	 *
-	 *                                  2* (1-n)/2n
-	 * mu2= -------------------------------------------------------------------
-	 *     2[ (du/dx)^2+(dv/dy)^2+1/4*(du/dy+dv/dx)^2+du/dx*dv/dy ]^[(3n-1)/2n]
-	 *
-	 *     where mu2 is the second viscosity, (u,v) the velocity 
-	 *     vector, and n the flow law exponent.
-	 *
-	 * If epsilon is NULL, it means this is the first time Gradjb is being run, and we 
-	 * return mu20, initial viscosity.
-	 */
-	
-	/*output: */
-	double viscosity_complement;
-
-	/*input strain rate: */
-	double exx,eyy,exy;
-
-	/*Intermediary value A and exponent e: */
-	double A,e;
-
-	if(epsilon){
-		exx=*(epsilon+0);
-		eyy=*(epsilon+1);
-		exy=*(epsilon+2);
-
-		/*Build viscosity: mu2=B/(2*A^e) */
-		A=pow(exx,2)+pow(eyy,2)+pow(exy,2)+exx*eyy;
-		if(A==0){
-			/*Maximum viscosity_complement for 0 shear areas: */
-			viscosity_complement=4.5*pow((double)10,(double)17);
-		}
-		else{
-			e=(n-1)/2/n;
-		
-			viscosity_complement=1/(2*pow(A,e));
-		}
-	}
-	else{
-		viscosity_complement=4.5*pow((double)10,(double)17);
-	}
-		
-	#ifdef _ISSM_DEBUG_
-	printf("viscosity_complement %lf\n",viscosity_complement);
-	#endif
-
-	/*Return: */
-	*pviscosity_complement=viscosity_complement;
-}
-
+/*}}}*/
+/*FUNCTION Matice::GetViscosity3d {{{1*/
 #undef __FUNCT__ 
 #define __FUNCT__ "MatIce::GetViscosity3d"
@@ -313,5 +271,6 @@
 	*pviscosity3d=viscosity3d;
 }
-
+/*}}}*/
+/*FUNCTION Matice::GetViscosity3dStokes {{{1*/
 #undef __FUNCT__ 
 #define __FUNCT__ "MatIce::GetViscosity3dStokes"
@@ -382,10 +341,81 @@
 	*pviscosity3d=viscosity3d;
 }
-double Matice::GetB(){
-	return B;
-}
-
-double Matice::GetN(){
-	return n;
-}
-
+/*}}}*/
+/*FUNCTION Matice::GetViscosityComplement {{{1*/
+#undef __FUNCT__ 
+#define __FUNCT__ "MatIce::GetViscosityComplement"
+void  Matice::GetViscosityComplement(double* pviscosity_complement, double* epsilon){
+
+	/*Return viscosity accounting for steady state power law creep [Thomas and MacAyeal, 1982]: 
+	 *
+	 *                                  2* (1-n)/2n
+	 * mu2= -------------------------------------------------------------------
+	 *     2[ (du/dx)^2+(dv/dy)^2+1/4*(du/dy+dv/dx)^2+du/dx*dv/dy ]^[(3n-1)/2n]
+	 *
+	 *     where mu2 is the second viscosity, (u,v) the velocity 
+	 *     vector, and n the flow law exponent.
+	 *
+	 * If epsilon is NULL, it means this is the first time Gradjb is being run, and we 
+	 * return mu20, initial viscosity.
+	 */
+	
+	/*output: */
+	double viscosity_complement;
+
+	/*input strain rate: */
+	double exx,eyy,exy;
+
+	/*Intermediary value A and exponent e: */
+	double A,e;
+
+	if(epsilon){
+		exx=*(epsilon+0);
+		eyy=*(epsilon+1);
+		exy=*(epsilon+2);
+
+		/*Build viscosity: mu2=B/(2*A^e) */
+		A=pow(exx,2)+pow(eyy,2)+pow(exy,2)+exx*eyy;
+		if(A==0){
+			/*Maximum viscosity_complement for 0 shear areas: */
+			viscosity_complement=4.5*pow((double)10,(double)17);
+		}
+		else{
+			e=(n-1)/2/n;
+		
+			viscosity_complement=1/(2*pow(A,e));
+		}
+	}
+	else{
+		viscosity_complement=4.5*pow((double)10,(double)17);
+	}
+		
+	#ifdef _ISSM_DEBUG_
+	printf("viscosity_complement %lf\n",viscosity_complement);
+	#endif
+
+	/*Return: */
+	*pviscosity_complement=viscosity_complement;
+}
+/*}}}*/
+/*FUNCTION Matice::MyRank {{{1*/
+int    Matice::MyRank(void){ 
+	extern int my_rank;
+	return my_rank; 
+}
+/*}}}*/
+/*FUNCTION Matice::SetB {{{1*/
+#undef __FUNCT__ 
+#define __FUNCT__ "Matice::SetB"
+void  Matice::SetB(double B_param){
+	B=B_param;
+}
+/*}}}*/
+/*FUNCTION Matice::UpdateFromInputs {{{1*/
+#undef __FUNCT__ 
+#define __FUNCT__ "Matice::UpdateFromInputs"
+void  Matice::UpdateFromInputs(void* inputs){
+	
+	//throw ErrorException(__FUNCT__," not supported yet!");
+
+}
+/*}}}*/
Index: /issm/trunk/src/c/objects/Matpar.cpp
===================================================================
--- /issm/trunk/src/c/objects/Matpar.cpp	(revision 2906)
+++ /issm/trunk/src/c/objects/Matpar.cpp	(revision 2907)
@@ -17,7 +17,11 @@
 
 		
+/*Object constructors and destructor*/
+/*FUNCTION Matpar::constructor {{{1*/
 Matpar::Matpar(){
 	return;
 }
+/*}}}1*/
+/*FUNCTION Matpar::creation {{{1*/
 Matpar::Matpar(int	matpar_mid, double	matpar_rho_ice, double	matpar_rho_water, double  matpar_heatcapacity, double  matpar_thermalconductivity, double  matpar_latentheat, double  matpar_beta, double  matpar_meltingpoint, double  matpar_mixed_layer_capacity, double  matpar_thermal_exchange_velocity, double  matpar_g){
 
@@ -37,9 +41,98 @@
 	return;
 }
-
+/*}}}1*/
+/*FUNCTION Matpar::destructor{{{1*/
 Matpar::~Matpar(){
 	return;
 }
-void Matpar::Echo(void){
+/*}}}1*/
+
+/*Object marshall*/
+/*FUNCTION Matpar::Demarshall {{{1*/
+void  Matpar::Demarshall(char** pmarshalled_dataset){
+
+	char* marshalled_dataset=NULL;
+
+	/*recover marshalled_dataset: */
+	marshalled_dataset=*pmarshalled_dataset;
+
+	/*this time, no need to get enum type, the pointer directly points to the beginning of the 
+	 *object data (thanks to DataSet::Demarshall):*/
+
+	memcpy(&mid,marshalled_dataset,sizeof(mid));marshalled_dataset+=sizeof(mid);
+	memcpy(&rho_ice,marshalled_dataset,sizeof(rho_ice));marshalled_dataset+=sizeof(rho_ice);
+	memcpy(&rho_water,marshalled_dataset,sizeof(rho_water));marshalled_dataset+=sizeof(rho_water);
+	memcpy(&heatcapacity,marshalled_dataset,sizeof(heatcapacity));marshalled_dataset+=sizeof(heatcapacity);
+	memcpy(&thermalconductivity,marshalled_dataset,sizeof(thermalconductivity));marshalled_dataset+=sizeof(thermalconductivity);
+	memcpy(&latentheat,marshalled_dataset,sizeof(latentheat));marshalled_dataset+=sizeof(latentheat);
+	memcpy(&beta,marshalled_dataset,sizeof(beta));marshalled_dataset+=sizeof(beta);
+	memcpy(&meltingpoint,marshalled_dataset,sizeof(meltingpoint));marshalled_dataset+=sizeof(meltingpoint);
+	memcpy(&mixed_layer_capacity,marshalled_dataset,sizeof(mixed_layer_capacity));marshalled_dataset+=sizeof(mixed_layer_capacity);
+	memcpy(&thermal_exchange_velocity,marshalled_dataset,sizeof(thermal_exchange_velocity));marshalled_dataset+=sizeof(thermal_exchange_velocity);
+	memcpy(&g,marshalled_dataset,sizeof(g));marshalled_dataset+=sizeof(g);
+
+	/*return: */
+	*pmarshalled_dataset=marshalled_dataset;
+	return;
+}
+/*}}}1*/
+/*FUNCTION Matpar::Marshall {{{1*/
+void  Matpar::Marshall(char** pmarshalled_dataset){
+
+	char* marshalled_dataset=NULL;
+	int   enum_type=0;
+
+	/*recover marshalled_dataset: */
+	marshalled_dataset=*pmarshalled_dataset;
+
+	/*get enum type of Matpar: */
+	enum_type=MatparEnum();
+	
+	/*marshall enum: */
+	memcpy(marshalled_dataset,&enum_type,sizeof(enum_type));marshalled_dataset+=sizeof(enum_type);
+	
+	/*marshall Matpar data: */
+	memcpy(marshalled_dataset,&mid,sizeof(mid));marshalled_dataset+=sizeof(mid);
+	memcpy(marshalled_dataset,&rho_ice,sizeof(rho_ice));marshalled_dataset+=sizeof(rho_ice);
+	memcpy(marshalled_dataset,&rho_water,sizeof(rho_water));marshalled_dataset+=sizeof(rho_water);
+	memcpy(marshalled_dataset,&heatcapacity,sizeof(heatcapacity));marshalled_dataset+=sizeof(heatcapacity);
+	memcpy(marshalled_dataset,&thermalconductivity,sizeof(thermalconductivity));marshalled_dataset+=sizeof(thermalconductivity);
+	memcpy(marshalled_dataset,&latentheat,sizeof(latentheat));marshalled_dataset+=sizeof(latentheat);
+	memcpy(marshalled_dataset,&beta,sizeof(beta));marshalled_dataset+=sizeof(beta);
+	memcpy(marshalled_dataset,&meltingpoint,sizeof(meltingpoint));marshalled_dataset+=sizeof(meltingpoint);
+	memcpy(marshalled_dataset,&mixed_layer_capacity,sizeof(mixed_layer_capacity));marshalled_dataset+=sizeof(mixed_layer_capacity);
+	memcpy(marshalled_dataset,&thermal_exchange_velocity,sizeof(thermal_exchange_velocity));marshalled_dataset+=sizeof(thermal_exchange_velocity);
+	memcpy(marshalled_dataset,&g,sizeof(g));marshalled_dataset+=sizeof(g);
+
+	*pmarshalled_dataset=marshalled_dataset;
+	return;
+}
+/*}}}1*/
+/*FUNCTION Matpar::MarshallSize {{{1*/
+int   Matpar::MarshallSize(){
+
+	return sizeof(mid)+
+		sizeof(rho_ice)+
+		sizeof(rho_water)+
+		sizeof(heatcapacity)+
+		sizeof(thermalconductivity)+
+		sizeof(latentheat)+
+		sizeof(beta)+
+		sizeof(meltingpoint)+
+		sizeof(mixed_layer_capacity)+
+		sizeof(thermal_exchange_velocity)+
+		sizeof(g)+
+		sizeof(int); //sizeof(int) for enum type
+}
+/*}}}1*/
+
+/*Object functions*/
+/*FUNCTION Matpar::copy {{{1*/
+Object* Matpar::copy() {
+	return new Matpar(*this); 
+}
+/*}}}1*/
+/*FUNCTION Matpar::DeepEcho {{{1*/
+void Matpar::DeepEcho(void){
 
 	printf("Matpar:\n");
@@ -56,7 +149,11 @@
 	printf("   g: %g\n",g);
 	return;
-}
-
-void Matpar::DeepEcho(void){
+}		
+/*}}}1*/
+/*FUNCTION Matpar::DistributeNumDofs {{{1*/
+void  Matpar::DistributeNumDofs(int* numdofspernode,int analysis_type){return;}
+/*}}}1*/
+/*FUNCTION Matpar::Echo {{{1*/
+void Matpar::Echo(void){
 
 	printf("Matpar:\n");
@@ -73,98 +170,79 @@
 	printf("   g: %g\n",g);
 	return;
-}		
-void  Matpar::Marshall(char** pmarshalled_dataset){
-
-	char* marshalled_dataset=NULL;
-	int   enum_type=0;
-
-	/*recover marshalled_dataset: */
-	marshalled_dataset=*pmarshalled_dataset;
-
-	/*get enum type of Matpar: */
-	enum_type=MatparEnum();
-	
-	/*marshall enum: */
-	memcpy(marshalled_dataset,&enum_type,sizeof(enum_type));marshalled_dataset+=sizeof(enum_type);
-	
-	/*marshall Matpar data: */
-	memcpy(marshalled_dataset,&mid,sizeof(mid));marshalled_dataset+=sizeof(mid);
-	memcpy(marshalled_dataset,&rho_ice,sizeof(rho_ice));marshalled_dataset+=sizeof(rho_ice);
-	memcpy(marshalled_dataset,&rho_water,sizeof(rho_water));marshalled_dataset+=sizeof(rho_water);
-	memcpy(marshalled_dataset,&heatcapacity,sizeof(heatcapacity));marshalled_dataset+=sizeof(heatcapacity);
-	memcpy(marshalled_dataset,&thermalconductivity,sizeof(thermalconductivity));marshalled_dataset+=sizeof(thermalconductivity);
-	memcpy(marshalled_dataset,&latentheat,sizeof(latentheat));marshalled_dataset+=sizeof(latentheat);
-	memcpy(marshalled_dataset,&beta,sizeof(beta));marshalled_dataset+=sizeof(beta);
-	memcpy(marshalled_dataset,&meltingpoint,sizeof(meltingpoint));marshalled_dataset+=sizeof(meltingpoint);
-	memcpy(marshalled_dataset,&mixed_layer_capacity,sizeof(mixed_layer_capacity));marshalled_dataset+=sizeof(mixed_layer_capacity);
-	memcpy(marshalled_dataset,&thermal_exchange_velocity,sizeof(thermal_exchange_velocity));marshalled_dataset+=sizeof(thermal_exchange_velocity);
-	memcpy(marshalled_dataset,&g,sizeof(g));marshalled_dataset+=sizeof(g);
-
-	*pmarshalled_dataset=marshalled_dataset;
-	return;
-}
-		
-int   Matpar::MarshallSize(){
-
-	return sizeof(mid)+
-		sizeof(rho_ice)+
-		sizeof(rho_water)+
-		sizeof(heatcapacity)+
-		sizeof(thermalconductivity)+
-		sizeof(latentheat)+
-		sizeof(beta)+
-		sizeof(meltingpoint)+
-		sizeof(mixed_layer_capacity)+
-		sizeof(thermal_exchange_velocity)+
-		sizeof(g)+
-		sizeof(int); //sizeof(int) for enum type
-}
-
+}
+/*}}}1*/
+/*FUNCTION Matpar::Enum {{{1*/
+int Matpar::Enum(void){
+
+	return MatparEnum();
+
+}
+/*}}}1*/
+/*FUNCTION Matpar::GetBeta {{{1*/
+double Matpar::GetBeta(){
+	return beta;
+}
+/*}}}1*/
+/*FUNCTION Matpar::GetG {{{1*/
+double Matpar::GetG(){
+	return g;
+}
+/*}}}1*/
+/*FUNCTION Matpar::GetHeatCapacity {{{1*/
+double Matpar::GetHeatCapacity(){
+	return heatcapacity;
+}
+/*}}}1*/
+/*FUNCTION Matpar::GetId {{{1*/
+int    Matpar::GetId(void){ return mid; }
+/*}}}1*/
+/*FUNCTION Matpar::GetLatentHeat {{{1*/
+double Matpar::GetLatentHeat(){
+	return latentheat;
+}
+/*}}}1*/
+/*FUNCTION Matpar::GetMeltingPoint {{{1*/
+double Matpar::GetMeltingPoint(){
+	return meltingpoint;
+}
+/*}}}1*/
+/*FUNCTION Matpar::GetMixedLayerCapacity {{{1*/
+double Matpar::GetMixedLayerCapacity(){
+	return mixed_layer_capacity;
+}
+/*}}}1*/
+/*FUNCTION Matpar::GetName {{{1*/
 char* Matpar::GetName(void){
 	return "matpar";
 }
-		
-
-void  Matpar::Demarshall(char** pmarshalled_dataset){
-
-	char* marshalled_dataset=NULL;
-
-	/*recover marshalled_dataset: */
-	marshalled_dataset=*pmarshalled_dataset;
-
-	/*this time, no need to get enum type, the pointer directly points to the beginning of the 
-	 *object data (thanks to DataSet::Demarshall):*/
-
-	memcpy(&mid,marshalled_dataset,sizeof(mid));marshalled_dataset+=sizeof(mid);
-	memcpy(&rho_ice,marshalled_dataset,sizeof(rho_ice));marshalled_dataset+=sizeof(rho_ice);
-	memcpy(&rho_water,marshalled_dataset,sizeof(rho_water));marshalled_dataset+=sizeof(rho_water);
-	memcpy(&heatcapacity,marshalled_dataset,sizeof(heatcapacity));marshalled_dataset+=sizeof(heatcapacity);
-	memcpy(&thermalconductivity,marshalled_dataset,sizeof(thermalconductivity));marshalled_dataset+=sizeof(thermalconductivity);
-	memcpy(&latentheat,marshalled_dataset,sizeof(latentheat));marshalled_dataset+=sizeof(latentheat);
-	memcpy(&beta,marshalled_dataset,sizeof(beta));marshalled_dataset+=sizeof(beta);
-	memcpy(&meltingpoint,marshalled_dataset,sizeof(meltingpoint));marshalled_dataset+=sizeof(meltingpoint);
-	memcpy(&mixed_layer_capacity,marshalled_dataset,sizeof(mixed_layer_capacity));marshalled_dataset+=sizeof(mixed_layer_capacity);
-	memcpy(&thermal_exchange_velocity,marshalled_dataset,sizeof(thermal_exchange_velocity));marshalled_dataset+=sizeof(thermal_exchange_velocity);
-	memcpy(&g,marshalled_dataset,sizeof(g));marshalled_dataset+=sizeof(g);
-
-	/*return: */
-	*pmarshalled_dataset=marshalled_dataset;
-	return;
-}
-
-int Matpar::Enum(void){
-
-	return MatparEnum();
-
-}
-
-int    Matpar::GetId(void){ return mid; }
-
+/*}}}1*/
+/*FUNCTION Matpar::GetRhoIce {{{1*/
+double Matpar::GetRhoIce(){
+	
+	return rho_ice;
+}
+/*}}}1*/
+/*FUNCTION Matpar::GetRhoWater {{{1*/
+double Matpar::GetRhoWater(){
+	return rho_water;
+}
+/*}}}1*/
+/*FUNCTION Matpar::GetThermalConductivity {{{1*/
+double Matpar::GetThermalConductivity(){
+	return thermalconductivity;
+}
+/*}}}1*/
+/*FUNCTION Matpar::GetThermalExchangeVelocity {{{1*/
+double Matpar::GetThermalExchangeVelocity(){
+	return thermal_exchange_velocity;
+}
+/*}}}1*/
+/*FUNCTION Matpar::MyRank {{{1*/
 int    Matpar::MyRank(void){ 
 	extern int my_rank;
 	return my_rank; 
 }
-void  Matpar::DistributeNumDofs(int* numdofspernode,int analysis_type){return;}
-
+/*}}}1*/
+/*FUNCTION Matpar::UpdateFromInputs {{{1*/
 #undef __FUNCT__ 
 #define __FUNCT__ "Matpar::UpdateFromInputs"
@@ -174,44 +252,3 @@
 
 }
-
-double Matpar::GetG(){
-	return g;
-}
-double Matpar::GetRhoIce(){
-	
-	return rho_ice;
-}
-double Matpar::GetRhoWater(){
-	return rho_water;
-}
-
-Object* Matpar::copy() {
-	return new Matpar(*this); 
-}
-
-double Matpar::GetMixedLayerCapacity(){
-	return mixed_layer_capacity;
-}
-
-double Matpar::GetThermalExchangeVelocity(){
-	return thermal_exchange_velocity;
-}
-
-double Matpar::GetHeatCapacity(){
-	return heatcapacity;
-}
-		
-double Matpar::GetThermalConductivity(){
-	return thermalconductivity;
-}
-		
-double Matpar::GetLatentHeat(){
-	return latentheat;
-}
-double Matpar::GetBeta(){
-	return beta;
-}
-double Matpar::GetMeltingPoint(){
-	return meltingpoint;
-}
-
+/*}}}1*/
Index: /issm/trunk/src/c/objects/Model.cpp
===================================================================
--- /issm/trunk/src/c/objects/Model.cpp	(revision 2906)
+++ /issm/trunk/src/c/objects/Model.cpp	(revision 2907)
@@ -18,4 +18,6 @@
 #include "stdio.h"
 
+/*Object constructors and destructor*/
+/*FUNCTION Model::constructor {{{1*/
 Model::Model(){
 
@@ -24,5 +26,6 @@
 
 }
-
+/*}}}1*/
+/*FUNCTION Model::destructor {{{1*/
 Model::~Model(){
 
@@ -30,23 +33,8 @@
 	/*do not delete active, it just points to one of the femmodels!*/
 }
-
-#undef __FUNCT__
-#define __FUNCT__ "Model::Echo"
-
-void Model::Echo(void){
-
-	printf("Models echo: \n");
-	printf("   number of fem models: %i\n",femmodels->Size());
-
-}
-#undef __FUNCT__
-#define __FUNCT__ "Model::DeepEcho"
-
-void Model::DeepEcho(void){
-
-	femmodels->Echo();
-
-}
-		
+/*}}}1*/
+
+/*Object functions*/
+/*FUNCTION Model::AddFormulation(ConstDataHandle MODEL, int analysis_type,int sub_analysis_type){{{1*/
 #undef __FUNCT__
 #define __FUNCT__ "Model::AddFormulation"
@@ -122,5 +110,6 @@
 
 }
-
+/*}}}1*/
+/*FUNCTION Model::AddFormulation(ConstDataHandle MODEL, int analysis_type) {{{1*/
 #undef __FUNCT__
 #define __FUNCT__ "Model::AddFormulation"
@@ -200,5 +189,27 @@
 
 }
-
+/*}}}1*/
+/*FUNCTION Model::DeepEcho {{{1*/
+#undef __FUNCT__
+#define __FUNCT__ "Model::DeepEcho"
+
+void Model::DeepEcho(void){
+
+	femmodels->Echo();
+
+}
+/*}}}1*/
+/*FUNCTION Model::Echo {{{1*/
+#undef __FUNCT__
+#define __FUNCT__ "Model::Echo"
+
+void Model::Echo(void){
+
+	printf("Models echo: \n");
+	printf("   number of fem models: %i\n",femmodels->Size());
+
+}
+/*}}}1*/
+/*FUNCTION Model::FindParam(int* pparameter,char* parametername {{{1*/
 #undef __FUNCT__
 #define __FUNCT__ "Model::FindParam"
@@ -219,4 +230,6 @@
 
 }
+/*}}}1*/
+/*FUNCTION Model::FindParam(double* pparameter,char* parametername ){{{1*/
 int   Model::FindParam(double* pparameter,char* parametername){
 	
@@ -235,4 +248,6 @@
 
 }
+/*}}}1*/
+/*FUNCTION Model::FindParam(double** pparameter,int* pM, int *pN,char* parametername) {{{1*/
 int   Model::FindParam(double** pparameter,int* pM, int *pN,char* parametername){
 	
@@ -251,4 +266,6 @@
 
 }
+/*}}}1*/
+/*FUNCTION Model::FindParam(char** pparameter,char* parametername) {{{1*/
 int   Model::FindParam(char** pparameter,char* parametername){
 	
@@ -266,5 +283,6 @@
 
 }
-
+/*}}}1*/
+/*FUNCTION Model::FindParam(int* pparameter,char* parametername,int analysis_type,int sub_analysis_type) {{{1*/
 int   Model::FindParam(int* pparameter,char* parametername,int analysis_type,int sub_analysis_type){
 	
@@ -279,5 +297,6 @@
 	femmodel->FindParam(pparameter,parametername);
 }
-
+/*}}}1*/
+/*FUNCTION Model::FindParam(double* pparameter,char* parametername,int analysis_type,int sub_analysis_type) {{{1*/
 int   Model::FindParam(double* pparameter,char* parametername,int analysis_type,int sub_analysis_type){
 	
@@ -292,6 +311,6 @@
 	femmodel->FindParam(pparameter,parametername);
 }
-
-
+/*}}}1*/
+/*FUNCTION Model::FindParam(double** pparameter,int* pM, int* pN,char* parametername,int analysis_type,int sub_analysis_type) {{{1*/
 int   Model::FindParam(double** pparameter,int* pM, int* pN,char* parametername,int analysis_type,int sub_analysis_type){
 
@@ -306,5 +325,6 @@
 	femmodel->FindParam(pparameter,pM, pN,parametername);
 }
-
+/*}}}1*/
+/*FUNCTION Model::FindParam(char** pparameter,char* parametername,int analysis_type,int sub_analysis_type) {{{1*/
 int   Model::FindParam(char** pparameter,char* parametername,int analysis_type,int sub_analysis_type){
 
@@ -319,5 +339,6 @@
 	femmodel->FindParam(pparameter,parametername);
 }
-
+/*}}}1*/
+/*FUNCTION Model::FindParam(int* pparameter,char* parametername,int analysis_type) {{{1*/
 int   Model::FindParam(int* pparameter,char* parametername,int analysis_type){
 	
@@ -332,5 +353,6 @@
 	femmodel->FindParam(pparameter,parametername);
 }
-
+/*}}}1*/
+/*FUNCTION Model::FindParam(double* pparameter,char* parametername,int analysis_type) {{{1*/
 int   Model::FindParam(double* pparameter,char* parametername,int analysis_type){
 	
@@ -345,6 +367,6 @@
 	femmodel->FindParam(pparameter,parametername);
 }
-
-
+/*}}}1*/
+/*FUNCTION Model::FindParam(double** pparameter,int* pM, int* pN,char* parametername,int analysis_type) {{{1*/
 int   Model::FindParam(double** pparameter,int* pM, int* pN,char* parametername,int analysis_type){
 
@@ -359,5 +381,6 @@
 	femmodel->FindParam(pparameter,pM,pN,parametername);
 }
-
+/*}}}1*/
+/*FUNCTION Model::FindParam(char** pparameter,char* parametername,int analysis_type) {{{1*/
 int   Model::FindParam(char** pparameter,char* parametername,int analysis_type){
 
@@ -372,5 +395,9 @@
 	femmodel->FindParam(pparameter,parametername);
 }	
-
+/*}}}1*/
+/*FUNCTION Model::GetActiveFormulation {{{1*/
+FemModel* Model::GetActiveFormulation(){return active;}
+/*}}}1*/
+/*FUNCTION Model::GetFormulation(int analysis_type,int sub_analysis_type) {{{1*/
 #undef __FUNCT__
 #define __FUNCT__ "Model::GetFormulation"
@@ -402,5 +429,6 @@
 	else return femmodel;
 }
-		
+/*}}}1*/
+/*FUNCTION Model::GetFormulation(int analysis_type) {{{1*/
 #undef __FUNCT__
 #define __FUNCT__ "Model::GetFormulation"
@@ -431,11 +459,8 @@
 	else return femmodel;
 }
-		
-
-FemModel* Model::GetActiveFormulation(){return active;}
-		
-
+/*}}}1*/
+/*FUNCTION Model::SetActiveFormulation {{{1*/
 void* Model::SetActiveFormulation(FemModel* femmodel){
 	active=femmodel;
 }
-
+/*}}}1*/
Index: /issm/trunk/src/c/objects/Node.cpp
===================================================================
--- /issm/trunk/src/c/objects/Node.cpp	(revision 2906)
+++ /issm/trunk/src/c/objects/Node.cpp	(revision 2907)
@@ -18,7 +18,11 @@
 #include "../include/typedefs.h"
 
+/*Object constructors and destructor*/
+/*FUNCTION Node constructor {{{1*/
 Node::Node(){
 	return;
 }
+/*}}}*/
+/*FUNCTION Node constructor {{{1*/
 Node::Node(int node_id,int node_partitionborder,int node_numdofs, double node_x[3],double node_sigma,int node_onbed,int node_onsurface,int node_upper_node_id,int node_onshelf,int node_onsheet){
 
@@ -52,82 +56,13 @@
 	return;
 }
-
+/*}}}*/
+/*FUNCTION Node destructeur{{{1*/
 Node::~Node(){
 	return;
 }
-		
-void Node::Echo(void){
-
-	int i;
-
-	printf("Node:\n");
-	printf("   id: %i\n",id);
-	printf("   partitionborder: %i\n",partitionborder);
-	printf("   clone: %i\n",clone);
-	printf("   numberofdofs: %i\n",numberofdofs);
-	printf("   x=[%g,%g,%g]\n",x[0],x[1],x[2]);
-	printf("   sigma=%g\n",sigma);
-	printf("   onbed: %i\n",onbed);
-	printf("   onsurface: %i\n",onsurface);
-	printf("   onshelf: %i\n",onshelf);
-	printf("   onsheet: %i\n",onsheet);
-	printf("   upper_node_id=%i\n",upper_node_id);
-	printf("   upper_node_offset=%i\n",upper_node_offset);
-	printf("   doflist:|");
-	for(i=0;i<numberofdofs;i++){
-		if(i>MAXDOFSPERNODE)break;
-		printf("%i|",doflist[i]);
-	}
-	printf("   doflist1:|");
-	printf("%i|",doflist1[0]);
-	printf("\n");
-
-	printf("   set membership: m,n,f,s sets \n");
-	for(i=0;i<numberofdofs;i++){
-		if(i>MAXDOFSPERNODE)break;
-		printf("      dof %i: %i %i %i %i\n",i,mset[i],nset[i],fset[i],sset[i]);
-	}
-	printf("\n");
-	if(upper_node)printf("   upper_node pointer: %p\n",upper_node);
-
-	return;
-}
-
-void Node::DeepEcho(void){
-
-	int i;
-
-	printf("Node:\n");
-	printf("   id: %i\n",id);
-	printf("   partitionborder: %i\n",partitionborder);
-	printf("   clone: %i\n",clone);
-	printf("   numberofdofs: %i\n",numberofdofs);
-	printf("   x=[%g,%g,%g]\n",x[0],x[1],x[2]);
-	printf("   sigma=%g\n",sigma);
-	printf("   onbed: %i\n",onbed);
-	printf("   onsurface: %i\n",onsurface);
-	printf("   onshelf: %i\n",onshelf);
-	printf("   onsheet: %i\n",onsheet);
-	printf("   upper_node_id=%i\n",upper_node_id);
-	printf("   upper_node_offset=%i\n",upper_node_offset);
-	printf("   doflist:|");
-	for(i=0;i<numberofdofs;i++){
-		if(i>MAXDOFSPERNODE)break;
-		printf("%i|",doflist[i]);
-	}
-	printf("   doflist1:|");
-	printf("%i|",doflist1[0]);
-	printf("\n");
-
-	printf("   set membership: m,n,f,s sets \n");
-	for(i=0;i<numberofdofs;i++){
-		if(i>MAXDOFSPERNODE)break;
-		printf("      dof %i: %i %i %i %i\n",i,mset[i],nset[i],fset[i],sset[i]);
-	}
-	printf("\n");
-	if(upper_node)printf("   upper_node pointer: %p\n",upper_node);
-
-	return;
-}		
+/*}}}*/
+
+/*Object marshall*/
+/*FUNCTION Node Marshall{{{1*/
 void  Node::Marshall(char** pmarshalled_dataset){
 
@@ -168,5 +103,6 @@
 	return;
 }
-		
+/*}}}*/
+/*FUNCTION Node MarshallSize{{{1*/
 int   Node::MarshallSize(){
 
@@ -192,10 +128,6 @@
 		sizeof(int); //sizeof(int) for enum type
 }
-
-char* Node::GetName(void){
-	return "node";
-}
-		
-
+/*}}}*/
+/*FUNCTION Node Demarshall{{{1*/
 void  Node::Demarshall(char** pmarshalled_dataset){
 
@@ -235,109 +167,51 @@
 	return;
 }
-
-int Node::Enum(void){
-
-	return NodeEnum();
-
-}
-int    Node::GetId(void){ return id; }
-
-int    Node::MyRank(void){ 
-	extern int my_rank;
-
-	return my_rank; 
-}
-
-void  Node::DistributeDofs(int* pdofcount,int* pdofcount1){
-
-	int i;
-	extern int my_rank;
-	int dofcount;
-	int dofcount1;
-
-	dofcount=*pdofcount;
-	dofcount1=*pdofcount1;
-	
-	if(clone){
-		/*This node is a clone! Don't distribute dofs, it will get them from another cpu!*/
-		return;
-	}
-
-	/*This node should distribute dofs, go ahead: */
-	for(i=0;i<numberofdofs;i++){
-		doflist[i]=dofcount+i;
-	}
-	dofcount+=numberofdofs;
-
-	doflist1[0]=dofcount1;
-	dofcount1+=1;
-
-	/*Assign output pointers: */
-	*pdofcount=dofcount;
-	*pdofcount1=dofcount1;
-
-	return; 
-}
-
-void  Node::UpdateDofs(int dofcount,int dofcount1){
-	
-	int i;
-	extern int my_rank;
-	
-	if(clone){
-		/*This node is a clone, don't update the dofs!: */
-		return;
-	}
-
-	/*This node should update the dofs, go ahead: */
-	for(i=0;i<numberofdofs;i++){
-		doflist[i]+=dofcount;
-	}
-	doflist1[0]+=dofcount1;
-
-	return; 
-}
-
-void  Node::ShowBorderDofs(int* borderdofs,int* borderdofs1){
-
-	int j;
-	extern int my_rank;
-	
-	/*Is this node on the partition border? */
-	if(!partitionborder)return;
-
-	/*Are we the cpu that created this node's dof list? */
-	if(clone)return;
-
-	/*Ok, we are on the partition border, and we did create the 
-	 * dofs for this node, plug the doflist into borderdofs: */
-	for(j=0;j<numberofdofs;j++){
-		*(borderdofs+numberofdofs*(id-1)+j)=doflist[j];
-	}
-	*(borderdofs1+(id-1)+0)=doflist1[0];
-
-	return;
-}
-
-void  Node::UpdateBorderDofs(int* allborderdofs,int* allborderdofs1){ 
-
-	int j;
-	extern int my_rank;
-	
-	/*Is this node on the partition border? */
-	if(!partitionborder)return;
-	
-	/*Are we the cpu that created this node's dof list? */
-	if(clone==0)return;
-
-	/*Ok, we are on the partition border, but we did not create 
-	 * the dofs for this node. Therefore, our doflist is garbage right 
-	 * now. Go pick it up in the allborderdofs: */
-	for(j=0;j<numberofdofs;j++){
-		doflist[j]=*(allborderdofs+numberofdofs*(id-1)+j);
-	}
-	doflist1[0]=*(allborderdofs1+(id-1)+0);
-	return; 
-}
+/*}}}*/
+		
+/*Object functions*/
+/*FUNCTION Node ApplyConstraints{{{1*/
+void  Node::ApplyConstraint(Vec yg,int dof,double value){
+
+	int index;
+
+	/*First, dof should be added in the s set, describing which 
+	 * dofs are constrained to a certain value (dirichlet boundary condition*/
+
+	DofInSSet(dof-1);
+
+	/*Second, we should add value into yg, at dof corresponding to doflist[dof], unless
+	 *  we are a clone!*/
+
+	if(!clone){
+
+		index=doflist[dof-1]; //matlab indexing
+
+		VecSetValues(yg,1,&index,&value,INSERT_VALUES);
+
+	}
+
+}
+/*}}}*/
+/*FUNCTION Node Configure {{{1*/
+#undef __FUNCT__ 
+#define __FUNCT__ "Node::Configure"
+void  Node::Configure(void* pnodes){
+
+	DataSet* nodes=NULL;
+
+	/*Recover pointers :*/
+	nodes=(DataSet*)pnodes;
+
+	/*Link this node with its upper node: */
+	ResolvePointers((Object**)&upper_node,&upper_node_id,&upper_node_offset,1,nodes);
+	
+}
+/*}}}*/
+/*FUNCTION Node copy {{{1*/
+Object* Node::copy() {
+	return new Node(*this); 
+}
+/*}}}*/
+/*FUNCTION Node CreatePartition{{{1*/
 void  Node::CreatePartition(Vec partition){ 
 
@@ -352,92 +226,6 @@
 	return;
 }
-		
-void  Node::SetClone(int* minranks){
-
-	extern int my_rank;
-
-	if (minranks[id-1]==my_rank){
-		clone=0;
-	}
-	else{
-		/*!there is a cpu with lower rank that has the same node, 
-		therefore, I am a clone*/
-		clone=1; 	
-	}
-
-}
-		
-
-int   Node::GetNumberOfDofs(){
-	
-	return numberofdofs;
-
-}
-
-int   Node::IsClone(){
-	
-	return clone;
-
-}
-
-
-void  Node::ApplyConstraint(Vec yg,int dof,double value){
-
-	int index;
-
-	/*First, dof should be added in the s set, describing which 
-	 * dofs are constrained to a certain value (dirichlet boundary condition*/
-
-	DofInSSet(dof-1);
-
-	/*Second, we should add value into yg, at dof corresponding to doflist[dof], unless
-	 *  we are a clone!*/
-
-	if(!clone){
-
-		index=doflist[dof-1]; //matlab indexing
-
-		VecSetValues(yg,1,&index,&value,INSERT_VALUES);
-
-	}
-
-}
-		
-
-void  Node::DofInSSet(int dof){
-
-	/*Put dof for this node into the s set (ie, this dof will be constrained 
-	 * to a fixed value during computations. */
-
-	mset[dof]=0; //m and n are mutually exclusive (m for rigid body modes)
-	nset[dof]=1;
-	fset[dof]=0; //n splits into f (for which we solve) and s (single point constraints)
-	sset[dof]=1;
-}
-
-void  Node::DofInMSet(int dof){
-
-	/*Put dof for this node into the m set (m set is for rigid body modes)*/
-
-	mset[dof]=1; //m and n are mutually exclusive (m for rigid body modes)
-	nset[dof]=0;
-	fset[dof]=0; //n splits into f (for which we solve) and s (single point constraints)
-	sset[dof]=0;
-}
-
-int  Node::DofIsInMSet(int dof){
-
-	if (mset[dof])return 1;
-	else return 0;
-
-}
-
-int   Node::GetDof(int dofindex){
-
-	return doflist[dofindex];
-
-}
-
-		
+/*}}}*/
+/*FUNCTION Node CreateVecSets {{{1*/
 void  Node::CreateVecSets(Vec pv_g,Vec pv_m,Vec pv_n,Vec pv_f,Vec pv_s){
 
@@ -472,88 +260,152 @@
 
 }
-		
-void  Node::GetDofList(int* outdoflist,int* pnumberofdofspernode){
+/*}}}*/
+/*FUNCTION Node DeepEcho{{{1*/
+void Node::DeepEcho(void){
 
 	int i;
+
+	printf("Node:\n");
+	printf("   id: %i\n",id);
+	printf("   partitionborder: %i\n",partitionborder);
+	printf("   clone: %i\n",clone);
+	printf("   numberofdofs: %i\n",numberofdofs);
+	printf("   x=[%g,%g,%g]\n",x[0],x[1],x[2]);
+	printf("   sigma=%g\n",sigma);
+	printf("   onbed: %i\n",onbed);
+	printf("   onsurface: %i\n",onsurface);
+	printf("   onshelf: %i\n",onshelf);
+	printf("   onsheet: %i\n",onsheet);
+	printf("   upper_node_id=%i\n",upper_node_id);
+	printf("   upper_node_offset=%i\n",upper_node_offset);
+	printf("   doflist:|");
 	for(i=0;i<numberofdofs;i++){
-		outdoflist[i]=doflist[i];
-	}
-	/*Assign output pointers:*/
-	*pnumberofdofspernode=numberofdofs;
-}
-
-int   Node::GetDofList1(void){
-	return doflist1[0];
-}
-
-double Node::GetX(){return x[0];}
-double Node::GetY(){return x[1];}
-double Node::GetZ(){return x[2];}
-double Node::GetSigma(){return sigma;}
-
-Object* Node::copy() {
-	return new Node(*this); 
-}
-
-#undef __FUNCT__ 
-#define __FUNCT__ "Node::UpdateFromInputs"
-void  Node::UpdateFromInputs(void* vinputs){
-	
-	ParameterInputs* inputs=NULL;
-	Node* node=this;
-	int dof[1]={0};
-
-	/*Recover parameter inputs: */
-	inputs=(ParameterInputs*)vinputs;
-
-	/*Update internal data if inputs holds new values: */
-	inputs->Recover("x",&x[0],1,dof,1,(void**)&node);
-	inputs->Recover("y",&x[1],1,dof,1,(void**)&node);
-	inputs->Recover("z",&x[2],1,dof,1,(void**)&node);
-
-	
-}
-
-#undef __FUNCT__ 
-#define __FUNCT__ "Node::Configure"
-void  Node::Configure(void* pnodes){
-
-	DataSet* nodes=NULL;
-
-	/*Recover pointers :*/
-	nodes=(DataSet*)pnodes;
-
-	/*Link this node with its upper node: */
-	ResolvePointers((Object**)&upper_node,&upper_node_id,&upper_node_offset,1,nodes);
-	
-}
-
-#undef __FUNCT__ 
-#define __FUNCT__ "Node::GetUpperNode"
-Node* Node::GetUpperNode(){
-	return upper_node;
-}
-		
-int   Node::IsOnBed(){
-	return onbed;
-}
-
-int   Node::IsOnSurface(){
-	return onsurface;
-}
-
-int   Node::IsOnShelf(){
-	return onshelf;
-}
-
-int   Node::IsOnSheet(){
-	return onsheet;
+		if(i>MAXDOFSPERNODE)break;
+		printf("%i|",doflist[i]);
+	}
+	printf("   doflist1:|");
+	printf("%i|",doflist1[0]);
+	printf("\n");
+
+	printf("   set membership: m,n,f,s sets \n");
+	for(i=0;i<numberofdofs;i++){
+		if(i>MAXDOFSPERNODE)break;
+		printf("      dof %i: %i %i %i %i\n",i,mset[i],nset[i],fset[i],sset[i]);
+	}
+	printf("\n");
+	if(upper_node)printf("   upper_node pointer: %p\n",upper_node);
+
+	return;
 }		
-void  Node::FreezeDof(int dof){
-	
-	DofInSSet(dof-1); //with 0 displacement for this dof.
-
-}
-		
+/*}}}*/
+/*FUNCTION Node DistributeDofs{{{1*/
+void  Node::DistributeDofs(int* pdofcount,int* pdofcount1){
+
+	int i;
+	extern int my_rank;
+	int dofcount;
+	int dofcount1;
+
+	dofcount=*pdofcount;
+	dofcount1=*pdofcount1;
+	
+	if(clone){
+		/*This node is a clone! Don't distribute dofs, it will get them from another cpu!*/
+		return;
+	}
+
+	/*This node should distribute dofs, go ahead: */
+	for(i=0;i<numberofdofs;i++){
+		doflist[i]=dofcount+i;
+	}
+	dofcount+=numberofdofs;
+
+	doflist1[0]=dofcount1;
+	dofcount1+=1;
+
+	/*Assign output pointers: */
+	*pdofcount=dofcount;
+	*pdofcount1=dofcount1;
+
+	return; 
+}
+/*}}}*/
+/*FUNCTION Node DofInMSet{{{1*/
+void  Node::DofInMSet(int dof){
+
+	/*Put dof for this node into the m set (m set is for rigid body modes)*/
+
+	mset[dof]=1; //m and n are mutually exclusive (m for rigid body modes)
+	nset[dof]=0;
+	fset[dof]=0; //n splits into f (for which we solve) and s (single point constraints)
+	sset[dof]=0;
+}
+/*}}}*/
+/*FUNCTION Node DofInSSet {{{1*/
+void  Node::DofInSSet(int dof){
+
+	/*Put dof for this node into the s set (ie, this dof will be constrained 
+	 * to a fixed value during computations. */
+
+	mset[dof]=0; //m and n are mutually exclusive (m for rigid body modes)
+	nset[dof]=1;
+	fset[dof]=0; //n splits into f (for which we solve) and s (single point constraints)
+	sset[dof]=1;
+}
+/*}}}*/
+/*FUNCTION Node DofIsInMSet{{{1*/
+int  Node::DofIsInMSet(int dof){
+
+	if (mset[dof])return 1;
+	else return 0;
+
+}
+/*}}}*/
+/*FUNCTION Node Echo{{{1*/
+void Node::Echo(void){
+
+	int i;
+
+	printf("Node:\n");
+	printf("   id: %i\n",id);
+	printf("   partitionborder: %i\n",partitionborder);
+	printf("   clone: %i\n",clone);
+	printf("   numberofdofs: %i\n",numberofdofs);
+	printf("   x=[%g,%g,%g]\n",x[0],x[1],x[2]);
+	printf("   sigma=%g\n",sigma);
+	printf("   onbed: %i\n",onbed);
+	printf("   onsurface: %i\n",onsurface);
+	printf("   onshelf: %i\n",onshelf);
+	printf("   onsheet: %i\n",onsheet);
+	printf("   upper_node_id=%i\n",upper_node_id);
+	printf("   upper_node_offset=%i\n",upper_node_offset);
+	printf("   doflist:|");
+	for(i=0;i<numberofdofs;i++){
+		if(i>MAXDOFSPERNODE)break;
+		printf("%i|",doflist[i]);
+	}
+	printf("   doflist1:|");
+	printf("%i|",doflist1[0]);
+	printf("\n");
+
+	printf("   set membership: m,n,f,s sets \n");
+	for(i=0;i<numberofdofs;i++){
+		if(i>MAXDOFSPERNODE)break;
+		printf("      dof %i: %i %i %i %i\n",i,mset[i],nset[i],fset[i],sset[i]);
+	}
+	printf("\n");
+	if(upper_node)printf("   upper_node pointer: %p\n",upper_node);
+
+	return;
+}
+/*}}}*/
+/*FUNCTION Node Enum{{{1*/
+int Node::Enum(void){
+
+	return NodeEnum();
+
+}
+/*}}}*/
+/*FUNCTION Node FieldDepthAverageAtBase{{{1*/
 #undef __FUNCT__ 
 #define __FUNCT__ "Node::FieldDepthAverageAtBase"
@@ -678,5 +530,6 @@
 	}
 }
-
+/*}}}*/
+/*FUNCTION Node FieldExtrude {{{1*/
 #undef __FUNCT__ 
 #define __FUNCT__ "Node::FieldExtrude"
@@ -789,5 +642,204 @@
 	} //if (extrude)
 }
-		
+/*}}}*/
+/*FUNCTION Node FreezeDof{{{1*/
+void  Node::FreezeDof(int dof){
+	
+	DofInSSet(dof-1); //with 0 displacement for this dof.
+
+}
+/*}}}*/
+/*FUNCTION Node GetDof {{{1*/
+int   Node::GetDof(int dofindex){
+
+	return doflist[dofindex];
+
+}
+/*}}}*/
+/*FUNCTION Node GetDofList{{{1*/
+void  Node::GetDofList(int* outdoflist,int* pnumberofdofspernode){
+
+	int i;
+	for(i=0;i<numberofdofs;i++){
+		outdoflist[i]=doflist[i];
+	}
+	/*Assign output pointers:*/
+	*pnumberofdofspernode=numberofdofs;
+}
+/*}}}*/
+/*FUNCTION Node GetDOfList1{{{1*/
+int   Node::GetDofList1(void){
+	return doflist1[0];
+}
+/*}}}*/
+/*FUNCTION Node GetId{{{1*/
+int    Node::GetId(void){ return id; }
+/*}}}*/
+/*FUNCTION Node GetName{{{1*/
+char* Node::GetName(void){
+	return "node";
+}
+/*}}}*/
+/*FUNCTION Node GetNumberOfDofs{{{1*/
+int   Node::GetNumberOfDofs(){
+	
+	return numberofdofs;
+
+}
+/*}}}*/
+/*FUNCTION Node GetSigma {{{1*/
+double Node::GetSigma(){return sigma;}
+/*}}}*/
+/*FUNCTION Node GetUpperNode {{{1*/
+#undef __FUNCT__ 
+#define __FUNCT__ "Node::GetUpperNode"
+Node* Node::GetUpperNode(){
+	return upper_node;
+}
+/*}}}*/
+/*FUNCTION Node GetX {{{1*/
+double Node::GetX(){return x[0];}
+/*}}}*/
+/*FUNCTION Node GetY {{{1*/
+double Node::GetY(){return x[1];}
+/*}}}*/
+/*FUNCTION Node GetZ {{{1*/
+double Node::GetZ(){return x[2];}
+/*}}}*/
+/*FUNCTION Node IsClone {{{1*/
+int   Node::IsClone(){
+	
+	return clone;
+
+}
+/*}}}*/
+/*FUNCTION Node IsOnBed {{{1*/
+int   Node::IsOnBed(){
+	return onbed;
+}
+/*}}}*/
+/*FUNCTION Node IsOnSheet {{{1*/
+int   Node::IsOnSheet(){
+	return onsheet;
+}		
+/*}}}*/
+/*FUNCTION Node IsOnShelf {{{1*/
+int   Node::IsOnShelf(){
+	return onshelf;
+}
+/*}}}*/
+/*FUNCTION Node IsOnSurface {{{1*/
+int   Node::IsOnSurface(){
+	return onsurface;
+}
+/*}}}*/
+/*FUNCTION Node MyRank{{{1*/
+int    Node::MyRank(void){ 
+	extern int my_rank;
+
+	return my_rank; 
+}
+/*}}}*/
+/*FUNCTION Node SetClone {{{1*/
+void  Node::SetClone(int* minranks){
+
+	extern int my_rank;
+
+	if (minranks[id-1]==my_rank){
+		clone=0;
+	}
+	else{
+		/*!there is a cpu with lower rank that has the same node, 
+		therefore, I am a clone*/
+		clone=1; 	
+	}
+
+}
+/*}}}*/
+/*FUNCTION Node ShowBorderDofs{{{1*/
+void  Node::ShowBorderDofs(int* borderdofs,int* borderdofs1){
+
+	int j;
+	extern int my_rank;
+	
+	/*Is this node on the partition border? */
+	if(!partitionborder)return;
+
+	/*Are we the cpu that created this node's dof list? */
+	if(clone)return;
+
+	/*Ok, we are on the partition border, and we did create the 
+	 * dofs for this node, plug the doflist into borderdofs: */
+	for(j=0;j<numberofdofs;j++){
+		*(borderdofs+numberofdofs*(id-1)+j)=doflist[j];
+	}
+	*(borderdofs1+(id-1)+0)=doflist1[0];
+
+	return;
+}
+/*}}}*/
+/*FUNCTION Node UpdateBorderDofs{{{1*/
+void  Node::UpdateBorderDofs(int* allborderdofs,int* allborderdofs1){ 
+
+	int j;
+	extern int my_rank;
+	
+	/*Is this node on the partition border? */
+	if(!partitionborder)return;
+	
+	/*Are we the cpu that created this node's dof list? */
+	if(clone==0)return;
+
+	/*Ok, we are on the partition border, but we did not create 
+	 * the dofs for this node. Therefore, our doflist is garbage right 
+	 * now. Go pick it up in the allborderdofs: */
+	for(j=0;j<numberofdofs;j++){
+		doflist[j]=*(allborderdofs+numberofdofs*(id-1)+j);
+	}
+	doflist1[0]=*(allborderdofs1+(id-1)+0);
+	return; 
+}
+/*}}}*/
+/*FUNCTION Node UpdateDofs{{{1*/
+void  Node::UpdateDofs(int dofcount,int dofcount1){
+	
+	int i;
+	extern int my_rank;
+	
+	if(clone){
+		/*This node is a clone, don't update the dofs!: */
+		return;
+	}
+
+	/*This node should update the dofs, go ahead: */
+	for(i=0;i<numberofdofs;i++){
+		doflist[i]+=dofcount;
+	}
+	doflist1[0]+=dofcount1;
+
+	return; 
+}
+/*}}}*/
+/*FUNCTION Node UpdateFromInputs {{{1*/
+#undef __FUNCT__ 
+#define __FUNCT__ "Node::UpdateFromInputs"
+void  Node::UpdateFromInputs(void* vinputs){
+	
+	ParameterInputs* inputs=NULL;
+	Node* node=this;
+	int dof[1]={0};
+
+	/*Recover parameter inputs: */
+	inputs=(ParameterInputs*)vinputs;
+
+	/*Update internal data if inputs holds new values: */
+	inputs->Recover("x",&x[0],1,dof,1,(void**)&node);
+	inputs->Recover("y",&x[1],1,dof,1,(void**)&node);
+	inputs->Recover("z",&x[2],1,dof,1,(void**)&node);
+
+	
+}
+/*}}}*/
+/*FUNCTION Node UpdateNodePosition {{{1*/
 void  Node::UpdateNodePosition(double* thickness,double* bed){
 
@@ -800,2 +852,3 @@
 
 }
+/*}}}*/
Index: /issm/trunk/src/c/objects/Numpar.cpp
===================================================================
--- /issm/trunk/src/c/objects/Numpar.cpp	(revision 2906)
+++ /issm/trunk/src/c/objects/Numpar.cpp	(revision 2907)
@@ -18,9 +18,11 @@
 
 
+/*Object constructors and destructor*/
+/*FUNCTION Numpar::constructor {{{1*/
 Numpar::Numpar(){
 	return;
 }
-
-
+/*}}}*/
+/*FUNCTION Numpar::creation {{{1*/
 Numpar::Numpar(int numpar_id){
 	id=numpar_id;
@@ -40,48 +42,13 @@
 	return;
 }
+/*}}}*/
+/*FUNCTION Numpar::destructor {{{1*/
 Numpar::~Numpar(){
 	return;
 }
-
-
-#undef __FUNCT__
-#define __FUNCT__ "Numpar::Echo"
-void Numpar::Echo(void){
-
-	printf("Numpar:\n");
-	printf("   id: %i\n",id);
-	printf("   meanvel: %g\n",meanvel);
-	printf("   epsvel: %g\n",epsvel);
-	printf("   artdiff: %i\n",artdiff);
-	printf("   viscosity_overshoot: %g\n",viscosity_overshoot);
-	printf("   stokesreconditioning: %g\n",stokesreconditioning);
-	printf("   control_type: %s\n",control_type);
-	printf("   cm_noisedmp: %g\n",cm_noisedmp);
-	printf("   cm_mindmp_value: %g\n",cm_mindmp_value);
-	printf("   cm_mindmp_slope: %g\n",cm_mindmp_slope);
-	printf("   cm_maxdmp_value: %g\n",cm_maxdmp_value);
-	printf("   cm_maxdmp_slope: %g\n",cm_maxdmp_slope);
-}
-
-#undef __FUNCT__
-#define __FUNCT__ "Numpar::DeepEcho"
-void Numpar::DeepEcho(void){
-
-	printf("Numpar:\n");
-	printf("   id: %i\n",id);
-	printf("   meanvel: %g\n",meanvel);
-	printf("   epsvel: %g\n",epsvel);
-	printf("   artdiff: %i\n",artdiff);
-	printf("   viscosity_overshoot: %g\n",viscosity_overshoot);
-	printf("   stokesreconditioning: %g\n",stokesreconditioning);
-	printf("   control_type: %s\n",control_type);
-	printf("   cm_noisedmp: %g\n",cm_noisedmp);
-	printf("   cm_mindmp_value: %g\n",cm_mindmp_value);
-	printf("   cm_mindmp_slope: %g\n",cm_mindmp_slope);
-	printf("   cm_maxdmp_value: %g\n",cm_maxdmp_value);
-	printf("   cm_maxdmp_slope: %g\n",cm_maxdmp_slope);
-}
-
-
+/*}}}*/
+
+/*Object marshall*/
+/*FUNCTION Numpar::Marshall {{{1*/
 void  Numpar::Marshall(char** pmarshalled_dataset){
 
@@ -115,5 +82,6 @@
 	return;
 }
-		
+/*}}}*/
+/*FUNCTION Numpar::MarshallSize{{{1*/
 int   Numpar::MarshallSize(){
 	return sizeof(id)
@@ -131,9 +99,6 @@
 		+sizeof(int); //sizeof(int) for enum type
 }
-
-char* Numpar::GetName(void){
-	return "beam";
-}
-
+/*}}}*/
+/*FUNCTION Numpar::Demarshall {{{1*/
 void  Numpar::Demarshall(char** pmarshalled_dataset){
 
@@ -164,17 +129,8 @@
 	return;
 }
-int Numpar::Enum(void){
-
-	return NumparEnum();
-
-}
-int    Numpar::GetId(void){ return id; }
-
-int    Numpar::MyRank(void){ 
-	extern int my_rank;
-	return my_rank; 
-}
-
-
+/*}}}*/
+
+/*Object functions*/
+/*FUNCTION Numpar::Configure {{{1*/
 #undef __FUNCT__ 
 #define __FUNCT__ "Numpar::Configure"
@@ -202,6 +158,74 @@
 	return;
 }
-
-
+/*}}}*/
+/*FUNCTION Numpar::copy {{{1*/
+Object* Numpar::copy() {
+	
+	return new Numpar(*this); 
+
+}
+/*}}}*/
+/*FUNCTION Numpar::DeepEcho {{{1*/
+#undef __FUNCT__
+#define __FUNCT__ "Numpar::DeepEcho"
+void Numpar::DeepEcho(void){
+
+	printf("Numpar:\n");
+	printf("   id: %i\n",id);
+	printf("   meanvel: %g\n",meanvel);
+	printf("   epsvel: %g\n",epsvel);
+	printf("   artdiff: %i\n",artdiff);
+	printf("   viscosity_overshoot: %g\n",viscosity_overshoot);
+	printf("   stokesreconditioning: %g\n",stokesreconditioning);
+	printf("   control_type: %s\n",control_type);
+	printf("   cm_noisedmp: %g\n",cm_noisedmp);
+	printf("   cm_mindmp_value: %g\n",cm_mindmp_value);
+	printf("   cm_mindmp_slope: %g\n",cm_mindmp_slope);
+	printf("   cm_maxdmp_value: %g\n",cm_maxdmp_value);
+	printf("   cm_maxdmp_slope: %g\n",cm_maxdmp_slope);
+}
+/*}}}*/
+/*FUNCTION Numpar::Echo{{{1*/
+#undef __FUNCT__
+#define __FUNCT__ "Numpar::Echo"
+void Numpar::Echo(void){
+
+	printf("Numpar:\n");
+	printf("   id: %i\n",id);
+	printf("   meanvel: %g\n",meanvel);
+	printf("   epsvel: %g\n",epsvel);
+	printf("   artdiff: %i\n",artdiff);
+	printf("   viscosity_overshoot: %g\n",viscosity_overshoot);
+	printf("   stokesreconditioning: %g\n",stokesreconditioning);
+	printf("   control_type: %s\n",control_type);
+	printf("   cm_noisedmp: %g\n",cm_noisedmp);
+	printf("   cm_mindmp_value: %g\n",cm_mindmp_value);
+	printf("   cm_mindmp_slope: %g\n",cm_mindmp_slope);
+	printf("   cm_maxdmp_value: %g\n",cm_maxdmp_value);
+	printf("   cm_maxdmp_slope: %g\n",cm_maxdmp_slope);
+}
+/*}}}*/
+/*FUNCTION Numpar::Enum {{{1*/
+int Numpar::Enum(void){
+
+	return NumparEnum();
+
+}
+/*}}}*/
+/*FUNCTION Numpar::GetId {{{1*/
+int    Numpar::GetId(void){ return id; }
+/*}}}*/
+/*FUNCTION Numpar::GetName {{{1*/
+char* Numpar::GetName(void){
+	return "beam";
+}
+/*}}}*/
+/*FUNCTION Numpar::MyRank {{{1*/
+int    Numpar::MyRank(void){ 
+	extern int my_rank;
+	return my_rank; 
+}
+/*}}}*/
+/*FUNCTION Numpar::UpdateFromInputs {{{1*/
 #undef __FUNCT__ 
 #define __FUNCT__ "Numpar::UpdateFromInputs"
@@ -227,9 +251,3 @@
 
 }
-		
-Object* Numpar::copy() {
-	
-	return new Numpar(*this); 
-
-}
-
+/*}}}*/
Index: /issm/trunk/src/c/objects/Param.cpp
===================================================================
--- /issm/trunk/src/c/objects/Param.cpp	(revision 2906)
+++ /issm/trunk/src/c/objects/Param.cpp	(revision 2907)
@@ -18,14 +18,6 @@
 
 		
-Param::Param(){
-	stringarray=NULL;
-	doublevec=NULL;
-	doublemat=NULL;
-	vec=NULL;
-	mat=NULL;
-
-	return;
-}
-
+/*Object constructors and destructor*/
+/*FUNCTION Param::constructor {{{1*/
 #undef __FUNCT__ 
 #define __FUNCT__  "Param constructor"
@@ -47,5 +39,17 @@
 
 }
-
+/*}}}*/
+/*FUNCTION Param::creation{{{1*/
+Param::Param(){
+	stringarray=NULL;
+	doublevec=NULL;
+	doublemat=NULL;
+	vec=NULL;
+	mat=NULL;
+
+	return;
+}
+/*}}}*/
+/*FUNCTION Param::destructor {{{1*/
 Param::~Param(){
 	
@@ -93,274 +97,8 @@
 	}
 }
-
-#undef __FUNCT__ 
-#define __FUNCT__  "Param echo"
-void Param::Echo(void){
-
-	int i,j;
-	
-	printf("Param:\n");
-	printf("   id: %i\n",id);
-	printf("   name: %s\n",name);
-	
-	switch(type){
-		case STRING:
-			printf("   string value: %s\n",string);
-			break;
-			
-		case  STRINGARRAY:
-			printf("   string array: %i strings\n",M);
-			for(i=0;i<M;i++){
-				printf("      %i: %s\n",i,stringarray[i]);
-			}
-			break;
-	
-		case DOUBLE:
-			printf("   double value: %g\n",ddouble);
-			break;
-		
-		case DOUBLEVEC:
-			/*printf("   double vector. size: %i ndof: %i\n",M,ndof);
-			for(i=0;i<M;i++)printf("%g\n",doublevec[i]);*/
-			break;
-	
-		case DOUBLEMAT:
-			printf("   double matrix. size: %i,%i\n",M,N);
-			for(i=0;i<M;i++){
-				for(j=0;j<N;j++){
-					printf("%g ",*(doublemat+N*i+j));
-				}
-				printf("\n");
-			}
-			break;
-
-		case PETSCVEC:
-			/*printf("   Petsc vector: \n");
-			VecView(vec,PETSC_VIEWER_STDOUT_WORLD);*/
-			break;
-
-		case  PETSCMAT:
-			/*printf("   Petsc matrix: \n");
-			MatView(mat,PETSC_VIEWER_STDOUT_WORLD);*/
-			break;
-
-		default:
-			throw ErrorException(__FUNCT__,exprintf("%s%i","unknow parameter type ",type));
-	}
-}
-
-#undef __FUNCT__ 
-#define __FUNCT__  "Param Deep Echo"
-void Param::DeepEcho(void){
-
-	int i,j;
-	
-	printf("Param:\n");
-	printf("   id: %i\n",id);
-	printf("   name: %s\n",name);
-	
-	switch(type){
-		case STRING:
-			printf("   string value: %s\n",string);
-			break;
-			
-		case  STRINGARRAY:
-			printf("   string array: %i strings\n",M);
-			for(i=0;i<M;i++){
-				printf("      %i: %s\n",i,stringarray[i]);
-			}
-	
-		case DOUBLE:
-			printf("   double value: %g\n",ddouble);
-			break;
-		
-		case DOUBLEVEC:
-			printf("   double vector. size: %i ndof: %i\n",M,ndof);
-			for(i=0;i<M;i++)printf("%g\n",doublevec[i]);
-			break;
-	
-		case DOUBLEMAT:
-			printf("   double matrix. size: %i,%i\n",M,N);
-			for(i=0;i<M;i++){
-				for(j=0;j<N;j++){
-					printf("%g ",*(doublemat+N*i+j));
-				}
-				printf("\n");
-			}
-			break;
-
-		case PETSCVEC:
-			printf("   Petsc vector: \n");
-			VecView(vec,PETSC_VIEWER_STDOUT_WORLD);
-			break;
-
-		case  PETSCMAT:
-			printf("   Petsc matrix: \n");
-			MatView(mat,PETSC_VIEWER_STDOUT_WORLD);
-			break;
-
-		default:
-			throw ErrorException(__FUNCT__,exprintf("%s%i","unknow parameter type ",type));
-	}
-}	
-#undef __FUNCT__
-#define __FUNCT__ "Marshall"
-
-void  Param::Marshall(char** pmarshalled_dataset){
-
-	char* marshalled_dataset=NULL;
-	int   enum_type=0;
-	double* serial_vec=NULL; 
-	double* serial_mat=NULL;
-	int i;
-	char* tempstring=NULL;
-
-	/*recover marshalled_dataset: */
-	marshalled_dataset=*pmarshalled_dataset;
-
-	/*get enum type of Param: */
-	enum_type=ParamEnum();
-	
-	/*marshall enum: */
-	memcpy(marshalled_dataset,&enum_type,sizeof(enum_type));marshalled_dataset+=sizeof(enum_type);
-	
-	/*marshall Param data: */
-	memcpy(marshalled_dataset,&id,sizeof(id));marshalled_dataset+=sizeof(id);
-	memcpy(marshalled_dataset,&name,sizeof(name));marshalled_dataset+=sizeof(name);
-	memcpy(marshalled_dataset,&type,sizeof(type));marshalled_dataset+=sizeof(type);
-	memcpy(marshalled_dataset,&ndof,sizeof(ndof));marshalled_dataset+=sizeof(ndof);
-
-	switch(type){
-		case STRING:
-			memcpy(marshalled_dataset,&string,sizeof(string));marshalled_dataset+=sizeof(string);
-			break;
-
-		case STRINGARRAY:
-			memcpy(marshalled_dataset,&M,sizeof(M));marshalled_dataset+=sizeof(M);
-			for(i=0;i<M;i++){
-				tempstring=stringarray[i];
-				int size=(strlen(tempstring)+1)*sizeof(char);
-				memcpy(marshalled_dataset,&size,sizeof(size));marshalled_dataset+=sizeof(size);
-				memcpy(marshalled_dataset,tempstring,size);marshalled_dataset+=size;
-			}
-			break;
-
-		case DOUBLE:
-			memcpy(marshalled_dataset,&ddouble,sizeof(ddouble));marshalled_dataset+=sizeof(ddouble);
-			break;
-		case DOUBLEVEC:
-			memcpy(marshalled_dataset,&M,sizeof(M));marshalled_dataset+=sizeof(M);
-			if(M){
-				memcpy(marshalled_dataset,doublevec,M*sizeof(double));
-				marshalled_dataset+=(M*sizeof(double));
-			}
-			break;
-
-		case DOUBLEMAT:
-			memcpy(marshalled_dataset,&M,sizeof(M));marshalled_dataset+=sizeof(M);
-			memcpy(marshalled_dataset,&N,sizeof(N));marshalled_dataset+=sizeof(N);
-			if(M*N){
-				memcpy(marshalled_dataset,doublemat,M*N*sizeof(double));
-				marshalled_dataset+=(M*N*sizeof(double));
-			}
-			break;
-
-		case PETSCVEC:
-			if(vec){
-				VecGetSize(vec,&M);
-				VecToMPISerial(&serial_vec,vec);
-				memcpy(marshalled_dataset,&M,sizeof(M));marshalled_dataset+=sizeof(M);
-				memcpy(marshalled_dataset,serial_vec,M*sizeof(double));marshalled_dataset+=(M*sizeof(double));
-				xfree((void**)&serial_vec);
-			}
-			else{
-				M=0;
-				memcpy(marshalled_dataset,&M,sizeof(M));marshalled_dataset+=sizeof(M);
-			}
-			break;
-
-		case PETSCMAT:
-			if(mat){
-				MatGetSize(mat,&M,&N);
-				MatToSerial(&serial_mat,mat);
-				memcpy(marshalled_dataset,&M,sizeof(M));marshalled_dataset+=sizeof(M);
-				memcpy(marshalled_dataset,&N,sizeof(N));marshalled_dataset+=sizeof(N);
-				memcpy(marshalled_dataset,serial_mat,M*N*sizeof(double));marshalled_dataset+=(M*N*sizeof(double));
-				xfree((void**)&serial_mat);
-			}
-			else{
-				M=0;
-				N=0;
-				memcpy(marshalled_dataset,&M,sizeof(M));marshalled_dataset+=sizeof(M);
-				memcpy(marshalled_dataset,&N,sizeof(N));marshalled_dataset+=sizeof(N);
-			}
-			break;
-		
-		default:
-			throw ErrorException(__FUNCT__,exprintf("%s%i","unknown parameter type",type));
-			break;
-	}
-
-	*pmarshalled_dataset=marshalled_dataset;
-	return;
-}
-		
-int   Param::MarshallSize(){
-
-	int size;
-	int i;
-	char* tempstring=NULL;
-
-	size=sizeof(id)+
-		sizeof(name)+
-		sizeof(type)+
-		sizeof(ndof)+
-		sizeof(int); //sizeof(int) for enum type
-
-	switch(type){
-		case STRING:
-			size+=sizeof(string);
-			break;
-
-		case STRINGARRAY:
-			size+=sizeof(M);
-			for(i=0;i<M;i++){
-				tempstring=stringarray[i];
-				size+=sizeof(M);
-				size+=(strlen(tempstring)+1)*sizeof(char);
-			}
-			break;
-
-		case DOUBLE:
-			size+= sizeof(ddouble);
-			break;
-		case DOUBLEVEC:
-			size+= sizeof(M)+M*sizeof(double);
-			break;
-
-		case DOUBLEMAT:
-			size+= sizeof(M)+sizeof(N)+M*N*sizeof(double);
-			break;
-
-		case PETSCVEC:
-			size+= sizeof(M)+M*sizeof(double);
-			break;
-
-		case PETSCMAT:
-			size+= sizeof(M)+sizeof(N)+M*N*sizeof(double);
-			break;
-
-		default:
-			throw ErrorException(__FUNCT__,exprintf("%s%i","unknown parameter type",type));
-	}
-
-	return size;
-}
-
-char* Param::GetName(void){
-	return "param";
-}
-		
-
+/*}}}*/
+
+/*Object marshall*/
+/*FUNCTION Param::Demarshall {{{1*/
 void  Param::Demarshall(char** pmarshalled_dataset){
 
@@ -481,15 +219,332 @@
 	return;
 }
-
-
+/*}}}*/
+/*FUNCTION Param::Marshall {{{1*/
+#undef __FUNCT__
+#define __FUNCT__ "Marshall"
+
+void  Param::Marshall(char** pmarshalled_dataset){
+
+	char* marshalled_dataset=NULL;
+	int   enum_type=0;
+	double* serial_vec=NULL; 
+	double* serial_mat=NULL;
+	int i;
+	char* tempstring=NULL;
+
+	/*recover marshalled_dataset: */
+	marshalled_dataset=*pmarshalled_dataset;
+
+	/*get enum type of Param: */
+	enum_type=ParamEnum();
+	
+	/*marshall enum: */
+	memcpy(marshalled_dataset,&enum_type,sizeof(enum_type));marshalled_dataset+=sizeof(enum_type);
+	
+	/*marshall Param data: */
+	memcpy(marshalled_dataset,&id,sizeof(id));marshalled_dataset+=sizeof(id);
+	memcpy(marshalled_dataset,&name,sizeof(name));marshalled_dataset+=sizeof(name);
+	memcpy(marshalled_dataset,&type,sizeof(type));marshalled_dataset+=sizeof(type);
+	memcpy(marshalled_dataset,&ndof,sizeof(ndof));marshalled_dataset+=sizeof(ndof);
+
+	switch(type){
+		case STRING:
+			memcpy(marshalled_dataset,&string,sizeof(string));marshalled_dataset+=sizeof(string);
+			break;
+
+		case STRINGARRAY:
+			memcpy(marshalled_dataset,&M,sizeof(M));marshalled_dataset+=sizeof(M);
+			for(i=0;i<M;i++){
+				tempstring=stringarray[i];
+				int size=(strlen(tempstring)+1)*sizeof(char);
+				memcpy(marshalled_dataset,&size,sizeof(size));marshalled_dataset+=sizeof(size);
+				memcpy(marshalled_dataset,tempstring,size);marshalled_dataset+=size;
+			}
+			break;
+
+		case DOUBLE:
+			memcpy(marshalled_dataset,&ddouble,sizeof(ddouble));marshalled_dataset+=sizeof(ddouble);
+			break;
+		case DOUBLEVEC:
+			memcpy(marshalled_dataset,&M,sizeof(M));marshalled_dataset+=sizeof(M);
+			if(M){
+				memcpy(marshalled_dataset,doublevec,M*sizeof(double));
+				marshalled_dataset+=(M*sizeof(double));
+			}
+			break;
+
+		case DOUBLEMAT:
+			memcpy(marshalled_dataset,&M,sizeof(M));marshalled_dataset+=sizeof(M);
+			memcpy(marshalled_dataset,&N,sizeof(N));marshalled_dataset+=sizeof(N);
+			if(M*N){
+				memcpy(marshalled_dataset,doublemat,M*N*sizeof(double));
+				marshalled_dataset+=(M*N*sizeof(double));
+			}
+			break;
+
+		case PETSCVEC:
+			if(vec){
+				VecGetSize(vec,&M);
+				VecToMPISerial(&serial_vec,vec);
+				memcpy(marshalled_dataset,&M,sizeof(M));marshalled_dataset+=sizeof(M);
+				memcpy(marshalled_dataset,serial_vec,M*sizeof(double));marshalled_dataset+=(M*sizeof(double));
+				xfree((void**)&serial_vec);
+			}
+			else{
+				M=0;
+				memcpy(marshalled_dataset,&M,sizeof(M));marshalled_dataset+=sizeof(M);
+			}
+			break;
+
+		case PETSCMAT:
+			if(mat){
+				MatGetSize(mat,&M,&N);
+				MatToSerial(&serial_mat,mat);
+				memcpy(marshalled_dataset,&M,sizeof(M));marshalled_dataset+=sizeof(M);
+				memcpy(marshalled_dataset,&N,sizeof(N));marshalled_dataset+=sizeof(N);
+				memcpy(marshalled_dataset,serial_mat,M*N*sizeof(double));marshalled_dataset+=(M*N*sizeof(double));
+				xfree((void**)&serial_mat);
+			}
+			else{
+				M=0;
+				N=0;
+				memcpy(marshalled_dataset,&M,sizeof(M));marshalled_dataset+=sizeof(M);
+				memcpy(marshalled_dataset,&N,sizeof(N));marshalled_dataset+=sizeof(N);
+			}
+			break;
+		
+		default:
+			throw ErrorException(__FUNCT__,exprintf("%s%i","unknown parameter type",type));
+			break;
+	}
+
+	*pmarshalled_dataset=marshalled_dataset;
+	return;
+}
+/*}}}*/
+/*FUNCTION Param::MarshallSize {{{1*/
+int   Param::MarshallSize(){
+
+	int size;
+	int i;
+	char* tempstring=NULL;
+
+	size=sizeof(id)+
+		sizeof(name)+
+		sizeof(type)+
+		sizeof(ndof)+
+		sizeof(int); //sizeof(int) for enum type
+
+	switch(type){
+		case STRING:
+			size+=sizeof(string);
+			break;
+
+		case STRINGARRAY:
+			size+=sizeof(M);
+			for(i=0;i<M;i++){
+				tempstring=stringarray[i];
+				size+=sizeof(M);
+				size+=(strlen(tempstring)+1)*sizeof(char);
+			}
+			break;
+
+		case DOUBLE:
+			size+= sizeof(ddouble);
+			break;
+		case DOUBLEVEC:
+			size+= sizeof(M)+M*sizeof(double);
+			break;
+
+		case DOUBLEMAT:
+			size+= sizeof(M)+sizeof(N)+M*N*sizeof(double);
+			break;
+
+		case PETSCVEC:
+			size+= sizeof(M)+M*sizeof(double);
+			break;
+
+		case PETSCMAT:
+			size+= sizeof(M)+sizeof(N)+M*N*sizeof(double);
+			break;
+
+		default:
+			throw ErrorException(__FUNCT__,exprintf("%s%i","unknown parameter type",type));
+	}
+
+	return size;
+}
+/*}}}*/
+
+/*Object functions*/
+/*FUNCTION Param::copy {{{1*/
+Object* Param::copy() {
+	return new Param(*this); 
+}
+/*}}}*/
+/*FUNCTION Param::Echo {{{1*/
+#undef __FUNCT__ 
+#define __FUNCT__  "Param echo"
+void Param::Echo(void){
+
+	int i,j;
+	
+	printf("Param:\n");
+	printf("   id: %i\n",id);
+	printf("   name: %s\n",name);
+	
+	switch(type){
+		case STRING:
+			printf("   string value: %s\n",string);
+			break;
+			
+		case  STRINGARRAY:
+			printf("   string array: %i strings\n",M);
+			for(i=0;i<M;i++){
+				printf("      %i: %s\n",i,stringarray[i]);
+			}
+			break;
+	
+		case DOUBLE:
+			printf("   double value: %g\n",ddouble);
+			break;
+		
+		case DOUBLEVEC:
+			/*printf("   double vector. size: %i ndof: %i\n",M,ndof);
+			for(i=0;i<M;i++)printf("%g\n",doublevec[i]);*/
+			break;
+	
+		case DOUBLEMAT:
+			printf("   double matrix. size: %i,%i\n",M,N);
+			for(i=0;i<M;i++){
+				for(j=0;j<N;j++){
+					printf("%g ",*(doublemat+N*i+j));
+				}
+				printf("\n");
+			}
+			break;
+
+		case PETSCVEC:
+			/*printf("   Petsc vector: \n");
+			VecView(vec,PETSC_VIEWER_STDOUT_WORLD);*/
+			break;
+
+		case  PETSCMAT:
+			/*printf("   Petsc matrix: \n");
+			MatView(mat,PETSC_VIEWER_STDOUT_WORLD);*/
+			break;
+
+		default:
+			throw ErrorException(__FUNCT__,exprintf("%s%i","unknow parameter type ",type));
+	}
+}
+/*}}}*/
+/*FUNCTION Param::Enum {{{1*/
 int Param::Enum(void){
 
 	return ParamEnum();
 }
-
+/*}}}*/
+/*FUNCTION Param::DeepEcho {{{1*/
+#undef __FUNCT__ 
+#define __FUNCT__  "Param Deep Echo"
+void Param::DeepEcho(void){
+
+	int i,j;
+	
+	printf("Param:\n");
+	printf("   id: %i\n",id);
+	printf("   name: %s\n",name);
+	
+	switch(type){
+		case STRING:
+			printf("   string value: %s\n",string);
+			break;
+			
+		case  STRINGARRAY:
+			printf("   string array: %i strings\n",M);
+			for(i=0;i<M;i++){
+				printf("      %i: %s\n",i,stringarray[i]);
+			}
+	
+		case DOUBLE:
+			printf("   double value: %g\n",ddouble);
+			break;
+		
+		case DOUBLEVEC:
+			printf("   double vector. size: %i ndof: %i\n",M,ndof);
+			for(i=0;i<M;i++)printf("%g\n",doublevec[i]);
+			break;
+	
+		case DOUBLEMAT:
+			printf("   double matrix. size: %i,%i\n",M,N);
+			for(i=0;i<M;i++){
+				for(j=0;j<N;j++){
+					printf("%g ",*(doublemat+N*i+j));
+				}
+				printf("\n");
+			}
+			break;
+
+		case PETSCVEC:
+			printf("   Petsc vector: \n");
+			VecView(vec,PETSC_VIEWER_STDOUT_WORLD);
+			break;
+
+		case  PETSCMAT:
+			printf("   Petsc matrix: \n");
+			MatView(mat,PETSC_VIEWER_STDOUT_WORLD);
+			break;
+
+		default:
+			throw ErrorException(__FUNCT__,exprintf("%s%i","unknow parameter type ",type));
+	}
+}	
+/*}}}*/
+/*FUNCTION Param::DistributeNumDofs {{{1*/
+void  Param::DistributeNumDofs(int* numdofspernode,int analysis_type){return;}
+/*}}}*/
+/*FUNCTION Param::GetM {{{1*/
+int   Param::GetM(){
+	return M;
+}
+/*}}}*/
+/*FUNCTION Param::GetN {{{1*/
+int   Param::GetN(){
+	return N;
+}
+/*}}}*/
+/*FUNCTION Param::GetNdof {{{1*/
+int   Param::GetNdof(){
+	return ndof;
+}
+/*}}}*/
+/*FUNCTION Param::GetParameterValue(Mat* pmat) {{{1*/
+#undef __FUNCT__ 
+#define __FUNCT__ "GetParameterValue"
+
+void  Param::GetParameterValue(Mat* pmat){
+	
+	Mat  outmat=NULL;
+	
+	if(type!=PETSCMAT)  throw ErrorException(__FUNCT__,exprintf("%s%i"," error message: trying to recover a Petsc mat from a param with type",type));
+	
+	if(mat){
+		MatDuplicate(mat,MAT_COPY_VALUES,&outmat);
+	}
+	*pmat=outmat;
+}
+/*}}}*/
+/*FUNCTION Param::GetName {{{1*/
+char* Param::GetName(void){
+	return "param";
+}
+/*}}}*/
+/*FUNCTION Param::GetParameterName {{{1*/
 char* Param::GetParameterName(void){ 
 	return name; 
 }
-
+/*}}}*/
+/*FUNCTION Param::GetParameterValue(double* pdouble) {{{1*/
 #undef __FUNCT__ 
 #define __FUNCT__ "GetParameterValue"
@@ -499,5 +554,6 @@
 	*pdouble=ddouble;
 }
-
+/*}}}*/
+/*FUNCTION Param::GetParameterValue(int* pinteger) {{{1*/
 #undef __FUNCT__ 
 #define __FUNCT__ "GetParameterValue"
@@ -507,5 +563,6 @@
 	*pinteger=(int)ddouble;
 }
-
+/*}}}*/
+/*FUNCTION Param::GetParameterValue(char** pstring) {{{1*/
 #undef __FUNCT__ 
 #define __FUNCT__ "GetParameterValue"
@@ -520,5 +577,6 @@
 
 }
-
+/*}}}*/
+/*FUNCTION Param::GetParameterValue(char*** pstringarray) {{{1*/
 #undef __FUNCT__ 
 #define __FUNCT__ "GetParameterValue"
@@ -537,5 +595,6 @@
 	*pstringarray=outstringarray;
 }
-
+/*}}}*/
+/*FUNCTION Param::GetParameterValue(double** pdoublearray) {{{1*/
 #undef __FUNCT__ 
 #define __FUNCT__ "GetParameterValue"
@@ -566,5 +625,6 @@
 	*pdoublearray=outdoublearray;
 }
-
+/*}}}*/
+/*FUNCTION Param::GetParameterValue(Vec* pvec) {{{1*/
 #undef __FUNCT__ 
 #define __FUNCT__ "GetParameterValue"
@@ -581,38 +641,20 @@
 	*pvec=outvec;
 }
-
-#undef __FUNCT__ 
-#define __FUNCT__ "GetParameterValue"
-
-void  Param::GetParameterValue(Mat* pmat){
-	
-	Mat  outmat=NULL;
-	
-	if(type!=PETSCMAT)  throw ErrorException(__FUNCT__,exprintf("%s%i"," error message: trying to recover a Petsc mat from a param with type",type));
-	
-	if(mat){
-		MatDuplicate(mat,MAT_COPY_VALUES,&outmat);
-	}
-	*pmat=outmat;
-}
-
+/*}}}*/
+/*FUNCTION Param::GetId {{{1*/
 int    Param::GetId(void){ return id; }
-
+/*}}}*/
+/*FUNCTION Param::GetType {{{1*/
+int   Param::GetType(){
+	return type;
+}
+/*}}}*/
+/*FUNCTION Param::MyRank {{{1*/
 int    Param::MyRank(void){ 
 	extern int my_rank;
 	return my_rank; 
 }
-
-void  Param::DistributeNumDofs(int* numdofspernode,int analysis_type){return;}
-		
-
-int   Param::GetType(){
-	return type;
-}
-
-Object* Param::copy() {
-	return new Param(*this); 
-}
-
+/*}}}*/
+/*FUNCTION Param::SetDouble(double value) {{{1*/
 #undef __FUNCT__
 #define __FUNCT__ "SetDouble"
@@ -621,5 +663,6 @@
 	ddouble=value;
 }
-
+/*}}}*/
+/*FUNCTION Param::SetDouble(int value) {{{1*/
 #undef __FUNCT__
 #define __FUNCT__ "SetDouble"
@@ -628,6 +671,50 @@
 	ddouble=(double)value;
 }
-
-
+/*}}}*/
+/*FUNCTION Param::SetDoubleMat {{{1*/
+#undef __FUNCT__
+#define __FUNCT__ "SetDoubleMat"
+void  Param::SetDoubleMat(double* value,int pM,int pN){
+	
+	if (type!=DOUBLEMAT) throw ErrorException(__FUNCT__,exprintf("%s%i"," trying to set doublematrix type",type));
+
+	this->M=pM;
+	this->N=pN;
+	if(this->M*this->N){
+		xfree((void**)&doublemat); doublemat=(double*)xcalloc(M*N,sizeof(double));
+		memcpy(doublemat,value,M*N*sizeof(double));
+	}
+}
+/*}}}*/
+/*FUNCTION Param::SetDoubleVec(double* value,int size) {{{1*/
+#undef __FUNCT__
+#define __FUNCT__ "SetDoubleVec"
+void  Param::SetDoubleVec(double* value,int size){
+	if (type!=DOUBLEVEC) throw ErrorException(__FUNCT__,exprintf("%s%i"," trying to set doublevecfor type",type));
+	
+	M=size;
+	if(M){
+		xfree((void**)&doublevec); doublevec=(double*)xmalloc(M*sizeof(double));
+		memcpy(doublevec,value,M*sizeof(double));
+	}
+	ndof=0; //this vector will not be repartitioned.
+
+}
+/*}}}*/
+/*FUNCTION Param::SetDoubleVec(double* value,int size, int numberofdofs) {{{1*/
+#undef __FUNCT__
+#define __FUNCT__ "SetDoubleVec"
+void  Param::SetDoubleVec(double* value,int size, int numberofdofs){
+	if (type!=DOUBLEVEC) throw ErrorException(__FUNCT__,exprintf("%s%i"," trying to set doublevecfor type",type));
+	
+	M=size;
+	if(M){
+		xfree((void**)&doublevec); doublevec=(double*)xmalloc(M*sizeof(double));
+		memcpy(doublevec,value,M*sizeof(double));
+	}
+	ndof=numberofdofs;
+}
+/*}}}*/
+/*FUNCTION Param::SetString {{{1*/
 #undef __FUNCT__
 #define __FUNCT__ "SetString"
@@ -636,5 +723,6 @@
 	strcpy(string,value);
 }
-
+/*}}}*/
+/*FUNCTION Param::SetStringArray {{{1*/
 #undef __FUNCT__
 #define __FUNCT__ "SetStringArray"
@@ -651,46 +739,6 @@
 	}
 }
-
-
-int   Param::GetM(){
-	return M;
-}
-
-int   Param::GetNdof(){
-	return ndof;
-}
-
-int   Param::GetN(){
-	return N;
-}
-
-#undef __FUNCT__
-#define __FUNCT__ "SetDoubleVec"
-void  Param::SetDoubleVec(double* value,int size){
-	if (type!=DOUBLEVEC) throw ErrorException(__FUNCT__,exprintf("%s%i"," trying to set doublevecfor type",type));
-	
-	M=size;
-	if(M){
-		xfree((void**)&doublevec); doublevec=(double*)xmalloc(M*sizeof(double));
-		memcpy(doublevec,value,M*sizeof(double));
-	}
-	ndof=0; //this vector will not be repartitioned.
-
-}
-
-#undef __FUNCT__
-#define __FUNCT__ "SetDoubleVec"
-void  Param::SetDoubleVec(double* value,int size, int numberofdofs){
-	if (type!=DOUBLEVEC) throw ErrorException(__FUNCT__,exprintf("%s%i"," trying to set doublevecfor type",type));
-	
-	M=size;
-	if(M){
-		xfree((void**)&doublevec); doublevec=(double*)xmalloc(M*sizeof(double));
-		memcpy(doublevec,value,M*sizeof(double));
-	}
-	ndof=numberofdofs;
-}
-
-
+/*}}}*/
+/*FUNCTION Param::SetVec {{{1*/
 #undef __FUNCT__
 #define __FUNCT__ "SetVec"
@@ -710,17 +758,3 @@
 
 }
-		
-
-#undef __FUNCT__
-#define __FUNCT__ "SetDoubleMat"
-void  Param::SetDoubleMat(double* value,int pM,int pN){
-	
-	if (type!=DOUBLEMAT) throw ErrorException(__FUNCT__,exprintf("%s%i"," trying to set doublematrix type",type));
-
-	this->M=pM;
-	this->N=pN;
-	if(this->M*this->N){
-		xfree((void**)&doublemat); doublemat=(double*)xcalloc(M*N,sizeof(double));
-		memcpy(doublemat,value,M*N*sizeof(double));
-	}
-}
+/*}}}*/
Index: /issm/trunk/src/c/objects/ParameterInputs.cpp
===================================================================
--- /issm/trunk/src/c/objects/ParameterInputs.cpp	(revision 2906)
+++ /issm/trunk/src/c/objects/ParameterInputs.cpp	(revision 2907)
@@ -18,4 +18,6 @@
 #include "./objects.h"
 
+/*Object constructors and destructor*/
+/*FUNCTION ParameterInputs::constructor {{{1*/
 ParameterInputs::ParameterInputs(){
 
@@ -23,5 +25,6 @@
 
 }
-		
+/*}}}*/
+/*FUNCTION ParameterInputs::destructor {{{1*/
 ParameterInputs::~ParameterInputs(){
 
@@ -29,5 +32,119 @@
 
 }
-
+/*}}}*/
+/*FUNCTION ParameterInputs::purge {{{1*/
+void ParameterInputs::purge(char* name){
+
+	int i;
+	Input* input=NULL;
+
+	/*Go through dataset, and figure out if an Input 
+	 * already has the name "name". If so, delete it : */
+	
+	for(i=0;i<dataset->Size();i++){
+		input=(Input*)dataset->GetObjectByOffset(i);
+
+		if (input->IsSameName(name)){
+			/*delete object: */
+			dataset->DeleteObject(input);
+		}
+	}
+}
+/*}}}*/
+
+/*Object functions*/
+/*FUNCTION ParameterInputs::Add(char* name,char* string) {{{1*/
+void ParameterInputs::Add(char* name,char* string){
+
+	Input* input=NULL;
+
+	/*First, purge object with same name: */
+	this->purge(name);
+
+	/*We are sure an input of the same name does not exist. Create new 
+	 * input: */
+	input=new Input(name,string);
+
+	/*Add input to dataset: */
+	dataset->AddObject(input);
+
+}
+/*}}}*/
+/*FUNCTION ParameterInputs::Add(char* name,int integer) {{{1*/
+void ParameterInputs::Add(char* name,int integer){
+	
+	Input* input=NULL;
+
+	/*First, purge object with same name: */
+	this->purge(name);
+
+	/*We are sure an input of the same name does not exist. Create new 
+	 * input: */
+	input=new Input(name,integer);
+
+	/*Add input to dataset: */
+	dataset->AddObject(input);
+}
+/*}}}*/
+/*FUNCTION ParameterInputs::Add(char* name,double scalar) {{{1*/
+void ParameterInputs::Add(char* name,double scalar){
+
+	Input* input=NULL;
+
+	/*First, purge object with same name: */
+	this->purge(name);
+
+	/*We are sure an input of the same name does not exist. Create new 
+	 * input: */
+	input=new Input(name,scalar);
+
+	/*Add input to dataset: */
+	dataset->AddObject(input);
+
+}
+/*}}}*/
+/*FUNCTION ParameterInputs::Add(char* name,double* vector,int ndof,int numberofnodes) {{{1*/
+void ParameterInputs::Add(char* name,double* vector,int ndof,int numberofnodes){
+	
+	Input* input=NULL;
+
+	/*First, purge object with same name: */
+	this->purge(name);
+
+	/*We are sure an input of the same name does not exist. Create new 
+	 * input: */
+	input=new Input(name,vector,ndof,numberofnodes);
+	
+	/*Add input to dataset: */
+	dataset->AddObject(input);
+
+}
+/*}}}*/
+/*FUNCTION ParameterInputs::Add(char* name,Vec petscvector,int ndof, int numberofnodes) {{{1*/
+void ParameterInputs::Add(char* name,Vec petscvector,int ndof, int numberofnodes){
+	
+	Input* input=NULL;
+
+	/*First, purge object with same name: */
+	this->purge(name);
+
+	/*We are sure an input of the same name does not exist. Create new 
+	 * input: */
+	input=new Input(name,petscvector,ndof,numberofnodes);
+
+	/*Add input to dataset: */
+	dataset->AddObject(input);
+
+}
+/*}}}*/
+/*FUNCTION ParameterInputs::DeepEcho {{{1*/
+void ParameterInputs::DeepEcho(){
+
+	printf("ParameterInputs echo: \n");
+	dataset->DeepEcho();
+
+}
+/*}}}*/
+/*FUNCTION ParameterInputs::Echo {{{1*/
 void ParameterInputs::Echo(){
 
@@ -36,109 +153,37 @@
 
 }
-
-void ParameterInputs::DeepEcho(){
-
-	printf("ParameterInputs echo: \n");
-	dataset->DeepEcho();
-
-}
-void ParameterInputs::purge(char* name){
-
-	int i;
-	Input* input=NULL;
-
-	/*Go through dataset, and figure out if an Input 
-	 * already has the name "name". If so, delete it : */
-	
-	for(i=0;i<dataset->Size();i++){
-		input=(Input*)dataset->GetObjectByOffset(i);
-
-		if (input->IsSameName(name)){
-			/*delete object: */
-			dataset->DeleteObject(input);
-		}
-	}
-}
-		
-
-void ParameterInputs::Add(char* name,char* string){
-
-	Input* input=NULL;
-
-	/*First, purge object with same name: */
-	this->purge(name);
-
-	/*We are sure an input of the same name does not exist. Create new 
-	 * input: */
-	input=new Input(name,string);
-
-	/*Add input to dataset: */
-	dataset->AddObject(input);
-
-}
-
-void ParameterInputs::Add(char* name,int integer){
-	
-	Input* input=NULL;
-
-	/*First, purge object with same name: */
-	this->purge(name);
-
-	/*We are sure an input of the same name does not exist. Create new 
-	 * input: */
-	input=new Input(name,integer);
-
-	/*Add input to dataset: */
-	dataset->AddObject(input);
-}
-
-
-void ParameterInputs::Add(char* name,double scalar){
-
-	Input* input=NULL;
-
-	/*First, purge object with same name: */
-	this->purge(name);
-
-	/*We are sure an input of the same name does not exist. Create new 
-	 * input: */
-	input=new Input(name,scalar);
-
-	/*Add input to dataset: */
-	dataset->AddObject(input);
-
-}
-
-void ParameterInputs::Add(char* name,double* vector,int ndof,int numberofnodes){
-	
-	Input* input=NULL;
-
-	/*First, purge object with same name: */
-	this->purge(name);
-
-	/*We are sure an input of the same name does not exist. Create new 
-	 * input: */
-	input=new Input(name,vector,ndof,numberofnodes);
-	
-	/*Add input to dataset: */
-	dataset->AddObject(input);
-
-}
-
-void ParameterInputs::Add(char* name,Vec petscvector,int ndof, int numberofnodes){
-	
-	Input* input=NULL;
-
-	/*First, purge object with same name: */
-	this->purge(name);
-
-	/*We are sure an input of the same name does not exist. Create new 
-	 * input: */
-	input=new Input(name,petscvector,ndof,numberofnodes);
-
-	/*Add input to dataset: */
-	dataset->AddObject(input);
-
-}
+/*}}}*/
+/*FUNCTION ParameterInputs::Get {{{1*/
+#undef __FUNCT__
+#define __FUNCT__ "ParameterInputs::Get"
+Vec ParameterInputs::Get(char* name,int* dofs, int numdofs){
+	
+	int i;
+	Input* input=NULL;
+	int found=0;
+	
+	/*output: */
+	Vec ug=NULL;
+
+	/*Go through dataset, and figure out if an Input 
+	 * has the name "name": */
+	for(i=0;i<dataset->Size();i++){
+		input=(Input*)dataset->GetObjectByOffset(i);
+
+		if (input->IsSameName(name)){
+			found=1;
+			break;
+		}
+	}
+
+	if(found==0)return NULL;
+	
+	/*call submethod: */
+	ug=input->Get(dofs,numdofs);
+
+	return ug;
+}
+/*}}}*/
+/*FUNCTION ParameterInputs::Recover(char* name, char** pstring) {{{1*/
 int ParameterInputs::Recover(char* name, char** pstring){
 
@@ -165,4 +210,6 @@
 
 }
+/*}}}*/
+/*FUNCTION ParameterInputs::Recover(char* name, int* pinteger) {{{1*/
 int ParameterInputs::Recover(char* name, int* pinteger){
 	
@@ -190,4 +237,6 @@
 	return found;
 }
+/*}}}*/
+/*FUNCTION ParameterInputs::Recover(char* name, double* pscalar) {{{1*/
 int ParameterInputs::Recover(char* name, double* pscalar){
 
@@ -216,5 +265,6 @@
 
 }
-
+/*}}}*/
+/*FUNCTION ParameterInputs::Recover(char* name,double* values, int ndof, int* dofs,int numnodes,void** pnodes) {{{1*/
 #undef __FUNCT__
 #define __FUNCT__ "ParameterInputs::Recover"
@@ -246,35 +296,6 @@
 	return 1;
 }
-		
-#undef __FUNCT__
-#define __FUNCT__ "ParameterInputs::Get"
-Vec ParameterInputs::Get(char* name,int* dofs, int numdofs){
-	
-	int i;
-	Input* input=NULL;
-	int found=0;
-	
-	/*output: */
-	Vec ug=NULL;
-
-	/*Go through dataset, and figure out if an Input 
-	 * has the name "name": */
-	for(i=0;i<dataset->Size();i++){
-		input=(Input*)dataset->GetObjectByOffset(i);
-
-		if (input->IsSameName(name)){
-			found=1;
-			break;
-		}
-	}
-
-	if(found==0)return NULL;
-	
-	/*call submethod: */
-	ug=input->Get(dofs,numdofs);
-
-	return ug;
-}
-
+/*}}}*/
+/*FUNCTION ParameterInputs::UpdateFromDakota {{{1*/
 #undef __FUNCT__
 #define __FUNCT__ "ParameterInputs::UpdateFromDakota"
@@ -359,2 +380,3 @@
 
 }
+/*}}}*/
Index: /issm/trunk/src/c/objects/Pengrid.cpp
===================================================================
--- /issm/trunk/src/c/objects/Pengrid.cpp	(revision 2906)
+++ /issm/trunk/src/c/objects/Pengrid.cpp	(revision 2907)
@@ -18,8 +18,11 @@
 
 		
+/*Object constructors and destructor*/
+/*FUNCTION Pengrid::constructor {{{1*/
 Pengrid::Pengrid(){
 	return;
 }
-
+/*}}}1*/
+/*FUNCTION Pengrid::creation {{{1*/
 Pengrid::Pengrid(int	pengrid_id, int pengrid_node_id,int pengrid_mparid, int pengrid_dof, int pengrid_active, double pengrid_penalty_offset,int pengrid_thermal_steadystate,int pengrid_stabilize_constraints){
 	
@@ -42,41 +45,13 @@
 	return;
 }
-
+/*}}}1*/
+/*FUNCTION Pengrid::destructor {{{1*/
 Pengrid::~Pengrid(){
 	return;
 }
-		
-void Pengrid::Echo(void){
-
-	printf("Pengrid:\n");
-	printf("   id: %i\n",id);
-	printf("   mparid: %i\n",mparid);
-	printf("   dof: %i\n",dof);
-	printf("   active: %i\n",active);
-	printf("   penalty_offset: %g\n",penalty_offset);
-	printf("   thermal_steadystate: %i\n",thermal_steadystate);
-	printf("   node_id: [%i]\n",node_id);
-	printf("   node_offset: [%i]\n",node_offset);
-	printf("   matpar_offset=%i\n",matpar_offset);
-	
-	return;
-}
-void Pengrid::DeepEcho(void){
-
-	printf("Pengrid:\n");
-	printf("   id: %i\n",id);
-	printf("   mparid: %i\n",mparid);
-	printf("   dof: %i\n",dof);
-	printf("   active: %i\n",active);
-	printf("   penalty_offset: %g\n",penalty_offset);
-	printf("   thermal_steadystate: %i\n",thermal_steadystate);
-	printf("   node_id: [%i]\n",node_id);
-	printf("   node_offset: [%i]\n",node_offset);
-	printf("   matpar_offset=%i\n",matpar_offset);
-	
-	if(node)node->Echo();
-	if(matpar)matpar->Echo();
-	return;
-}		
+/*}}}1*/
+		
+/*Object marshall*/
+/*FUNCTION Pengrid::Marshall {{{1*/
 void  Pengrid::Marshall(char** pmarshalled_dataset){
 
@@ -110,5 +85,6 @@
 	return;
 }
-		
+/*}}}1*/
+/*FUNCTION Pengrid::MarshallSize {{{1*/
 int   Pengrid::MarshallSize(){
 
@@ -127,10 +103,6 @@
 		sizeof(int); //sizeof(int) for enum type
 }
-
-char* Pengrid::GetName(void){
-	return "pengrid";
-}
-		
-
+/*}}}1*/
+/*FUNCTION Pengrid::Demarshall {{{1*/
 void  Pengrid::Demarshall(char** pmarshalled_dataset){
 
@@ -163,59 +135,238 @@
 	return;
 }
+/*}}}1*/
+
+/*Object functions*/
+/*FUNCTION Pengrid::copy {{{1*/
+Object* Pengrid::copy() {
+	return new Pengrid(*this); 
+}
+/*}}}1*/
+/*FUNCTION Pengrid::Configure {{{1*/
+#undef __FUNCT__ 
+#define __FUNCT__ "Pengrid::Configure"
+
+void  Pengrid::Configure(void* pelementsin,void* pnodesin,void* pmaterialsin){
+
+	DataSet* nodesin=NULL;
+	DataSet* materialsin=NULL;
+
+	/*Recover pointers :*/
+	nodesin=(DataSet*)pnodesin;
+	materialsin=(DataSet*)pmaterialsin;
+
+	/*Link this load with its nodes: */
+	ResolvePointers((Object**)&node,&node_id,&node_offset,1,nodesin);
+	ResolvePointers((Object**)&matpar,&mparid,&matpar_offset,1,materialsin);
+}
+/*}}}1*/
+/*FUNCTION Pengrid::CreateKMatrix {{{1*/
+#undef __FUNCT__ 
+#define __FUNCT__ "Pengrid::CreateKMatrix"
+
+void  Pengrid::CreateKMatrix(Mat Kgg,void* inputs,int analysis_type,int sub_analysis_type){
+
+	/*No loads applied, do nothing: */
+	return;
+
+}
+/*}}}1*/
+/*FUNCTION Pengrid::CreatePVector {{{1*/
+#undef __FUNCT__ 
+#define __FUNCT__ "Pengrid::CreatePVector"
+void  Pengrid::CreatePVector(Vec pg, void* inputs, int analysis_type,int sub_analysis_type){
+
+	/*No loads applied, do nothing: */
+	return;
+
+}
+/*}}}1*/
+/*FUNCTION Pengrid::DeepEcho {{{1*/
+void Pengrid::DeepEcho(void){
+
+	printf("Pengrid:\n");
+	printf("   id: %i\n",id);
+	printf("   mparid: %i\n",mparid);
+	printf("   dof: %i\n",dof);
+	printf("   active: %i\n",active);
+	printf("   penalty_offset: %g\n",penalty_offset);
+	printf("   thermal_steadystate: %i\n",thermal_steadystate);
+	printf("   node_id: [%i]\n",node_id);
+	printf("   node_offset: [%i]\n",node_offset);
+	printf("   matpar_offset=%i\n",matpar_offset);
+	
+	if(node)node->Echo();
+	if(matpar)matpar->Echo();
+	return;
+}		
+/*}}}1*/
+/*FUNCTION Pengrid::DistributenumDofs {{{1*/
+void  Pengrid::DistributeNumDofs(int* numdofpernode,int analysis_type,int sub_analysis_type){return;}
+/*}}}1*/
+/*FUNCTION Pengrid::Echo {{{1*/
+void Pengrid::Echo(void){
+
+	printf("Pengrid:\n");
+	printf("   id: %i\n",id);
+	printf("   mparid: %i\n",mparid);
+	printf("   dof: %i\n",dof);
+	printf("   active: %i\n",active);
+	printf("   penalty_offset: %g\n",penalty_offset);
+	printf("   thermal_steadystate: %i\n",thermal_steadystate);
+	printf("   node_id: [%i]\n",node_id);
+	printf("   node_offset: [%i]\n",node_offset);
+	printf("   matpar_offset=%i\n",matpar_offset);
+	
+	return;
+}
+/*}}}1*/
+/*FUNCTION Pengrid::Enum {{{1*/
 int Pengrid::Enum(void){
 
 	return PengridEnum();
 }
-
+/*}}}1*/
+/*FUNCTION Pengrid::GetDofList {{{1*/
+void  Pengrid::GetDofList(int* doflist,int* pnumberofdofspernode){
+
+	int j;
+	int doflist_per_node[MAXDOFSPERNODE];
+	int numberofdofspernode;
+	
+	node->GetDofList(&doflist_per_node[0],&numberofdofspernode);
+	for(j=0;j<numberofdofspernode;j++){
+		doflist[j]=doflist_per_node[j];
+	}
+
+	/*Assign output pointers:*/
+	*pnumberofdofspernode=numberofdofspernode;
+}
+/*}}}1*/
+/*FUNCTION Pengrid::GetId {{{1*/
 int    Pengrid::GetId(void){ return id; }
-
+/*}}}1*/
+/*FUNCTION Pengrid::GetName {{{1*/
+char* Pengrid::GetName(void){
+	return "pengrid";
+}
+/*}}}1*/
+/*FUNCTION Pengrid::MyRank {{{1*/
 int    Pengrid::MyRank(void){ 
 	extern int my_rank;
 	return my_rank; 
 }
-void  Pengrid::DistributeNumDofs(int* numdofpernode,int analysis_type,int sub_analysis_type){return;}
-
-#undef __FUNCT__ 
-#define __FUNCT__ "Pengrid::Configure"
-
-void  Pengrid::Configure(void* pelementsin,void* pnodesin,void* pmaterialsin){
-
-	DataSet* nodesin=NULL;
-	DataSet* materialsin=NULL;
-
-	/*Recover pointers :*/
-	nodesin=(DataSet*)pnodesin;
-	materialsin=(DataSet*)pmaterialsin;
-
-	/*Link this load with its nodes: */
-	ResolvePointers((Object**)&node,&node_id,&node_offset,1,nodesin);
-	ResolvePointers((Object**)&matpar,&mparid,&matpar_offset,1,materialsin);
-}
-
-
-#undef __FUNCT__ 
-#define __FUNCT__ "Pengrid::CreateKMatrix"
-
-void  Pengrid::CreateKMatrix(Mat Kgg,void* inputs,int analysis_type,int sub_analysis_type){
-
-	/*No loads applied, do nothing: */
-	return;
-
-}
-
-#undef __FUNCT__ 
-#define __FUNCT__ "Pengrid::CreatePVector"
-void  Pengrid::CreatePVector(Vec pg, void* inputs, int analysis_type,int sub_analysis_type){
-
-	/*No loads applied, do nothing: */
-	return;
-
-}
-#undef __FUNCT__ 
-#define __FUNCT__ "Pengrid::UpdateFromInputs"
-void  Pengrid::UpdateFromInputs(void* inputs){
-	
-}
-
+/*}}}1*/
+/*FUNCTION Pengrid::PenaltyConstrain {{{1*/
+#undef __FUNCT__ 
+#define __FUNCT__ "Pengrid::PenaltyConstrain"
+void  Pengrid::PenaltyConstrain(int* punstable,void* vinputs,int analysis_type,int sub_analysis_type){
+
+	if ((analysis_type==DiagnosticAnalysisEnum()) && ((sub_analysis_type==StokesAnalysisEnum()))){
+
+		/*No penalty to check*/
+		return;
+
+	}
+	else if (analysis_type==ThermalAnalysisEnum()){
+		
+		PenaltyConstrainThermal(punstable,vinputs,analysis_type,sub_analysis_type);
+		
+	}
+	else if (analysis_type==MeltingAnalysisEnum()){
+			
+		/*No penalty to check*/
+		return;
+
+	}
+	else{
+		throw ErrorException(__FUNCT__,exprintf("%s%i%s%i%s","analysis: ",analysis_type," and sub_analysis_type: ",sub_analysis_type," not supported yet"));
+	}
+
+}
+/*}}}1*/
+/*FUNCTION Pengrid::PenaltyConstrainThermal {{{1*/
+#undef __FUNCT__ 
+#define __FUNCT__ "Pengrid::PenaltyConstrainThermal"
+void  Pengrid::PenaltyConstrainThermal(int* punstable,void* vinputs,int analysis_type,int sub_analysis_type){
+
+	//   The penalty is stable if it doesn't change during to successive iterations.   
+
+	int found=0;
+	const int numgrids=1;
+
+
+	double pressure;
+	double temperature;
+	double beta,t_pmp;
+	double meltingpoint;
+	int    new_active;
+	int  dofs1[1]={0};
+	int  unstable=0;
+	int  reset_penalties=0;
+	
+	ParameterInputs* inputs=NULL;
+
+	/*check that pengrid is not a clone (penalty to be added only once)*/
+	if (node->IsClone()){
+		unstable=0;
+		*punstable=unstable;
+		return;
+	}
+
+	/*recover pointers: */
+	inputs=(ParameterInputs*)vinputs;
+
+	//First recover beta, pressure and temperature vectors;
+	found=inputs->Recover("pressure",&pressure,1,dofs1,numgrids,(void**)&node);
+	if(!found)throw ErrorException(__FUNCT__," could not find pressure in inputs!");
+
+	found=inputs->Recover("temperature",&temperature,1,dofs1,numgrids,(void**)&node);
+	if(!found)throw ErrorException(__FUNCT__," could not find temperature in inputs!");
+	
+	found=inputs->Recover("reset_penalties",&reset_penalties);
+
+	if(reset_penalties)zigzag_counter=0;
+
+	//Compute pressure melting point
+	meltingpoint=matpar->GetMeltingPoint();
+	beta=matpar->GetBeta();
+
+	t_pmp=meltingpoint-beta*pressure;
+
+	//Figure out if temperature is over melting_point, in which case, this penalty needs to be activated.
+
+	if (temperature>t_pmp){
+		new_active=1;
+	}
+	else{
+		new_active=0;
+	}
+
+
+	//Figure out stability of this penalty
+	if (active==new_active){
+		unstable=0;
+	}
+	else{
+		unstable=1;
+		if(stabilize_constraints)zigzag_counter++;
+	}
+
+	/*If penalty keeps zigzagging more than 5 times: */
+	if(stabilize_constraints){
+		if(zigzag_counter>stabilize_constraints){
+			unstable=0;
+			active=1;
+		}
+	}
+
+	//Set penalty flag
+	active=new_active;
+
+	//*Assign output pointers:*/
+	*punstable=unstable;
+}
+/*}}}1*/
+/*FUNCTION Pengrid::PenaltyCreateMatrix {{{1*/
 #undef __FUNCT__ 
 #define __FUNCT__ "Pengrid::PenaltyCreateKMatrix"
@@ -241,5 +392,6 @@
 
 }
-
+/*}}}1*/
+/*FUNCTION Pengrid::PenaltyCreateKMatrixDiagnosticStokes {{{1*/
 #undef __FUNCT__ 
 #define __FUNCT__ "Pengrid::PenaltyCreateKMatrixDiagnosticStokes"
@@ -280,35 +432,6 @@
 	MatSetValues(Kgg,numdof,doflist,numdof,doflist,(const double*)Ke,ADD_VALUES);
 }
-
-#undef __FUNCT__ 
-#define __FUNCT__ "Pengrid::PenaltyCreateKMatrixThermal"
-void  Pengrid::PenaltyCreateKMatrixThermal(Mat Kgg,void* vinputs,double kmax,int analysis_type,int sub_analysis_type){
-
-	int found=0;
-	
-	const int numgrids=1;
-	const int NDOF1=1;
-	const int numdof=numgrids*NDOF1;
-	double Ke[numdof][numdof];
-	int       doflist[numdof];
-	int       numberofdofspernode;
-
-	ParameterInputs* inputs=NULL;
-
-	/*recover pointers: */
-	inputs=(ParameterInputs*)vinputs;
-
-
-	if(!active)return;
-
-	/*Get dof list: */
-	GetDofList(&doflist[0],&numberofdofspernode);
-
-	Ke[0][0]=kmax*pow((double)10,penalty_offset);
-	
-	/*Add Ke to global matrix Kgg: */
-	MatSetValues(Kgg,numdof,doflist,numdof,doflist,(const double*)Ke,ADD_VALUES);
-}
-
+/*}}}1*/
+/*FUNCTION Pengrid::PenaltyCreateKMatrixMelting {{{1*/
 #undef __FUNCT__ 
 #define __FUNCT__ "Pengrid::PenaltyCreateKMatrixMelting"
@@ -359,5 +482,37 @@
 	MatSetValues(Kgg,numdof,doflist,numdof,doflist,(const double*)Ke,ADD_VALUES);
 }
-
+/*}}}1*/
+/*FUNCTION Pengrid::PenaltyCreateKMatrixThermal {{{1*/
+#undef __FUNCT__ 
+#define __FUNCT__ "Pengrid::PenaltyCreateKMatrixThermal"
+void  Pengrid::PenaltyCreateKMatrixThermal(Mat Kgg,void* vinputs,double kmax,int analysis_type,int sub_analysis_type){
+
+	int found=0;
+	
+	const int numgrids=1;
+	const int NDOF1=1;
+	const int numdof=numgrids*NDOF1;
+	double Ke[numdof][numdof];
+	int       doflist[numdof];
+	int       numberofdofspernode;
+
+	ParameterInputs* inputs=NULL;
+
+	/*recover pointers: */
+	inputs=(ParameterInputs*)vinputs;
+
+
+	if(!active)return;
+
+	/*Get dof list: */
+	GetDofList(&doflist[0],&numberofdofspernode);
+
+	Ke[0][0]=kmax*pow((double)10,penalty_offset);
+	
+	/*Add Ke to global matrix Kgg: */
+	MatSetValues(Kgg,numdof,doflist,numdof,doflist,(const double*)Ke,ADD_VALUES);
+}
+/*}}}1*/
+/*FUNCTION Pengrid::PenaltyCreatePVector {{{1*/
 #undef __FUNCT__ 
 #define __FUNCT__ "Pengrid::PenaltyCreatePVector"
@@ -385,66 +540,6 @@
 
 }
-
-Object* Pengrid::copy() {
-	return new Pengrid(*this); 
-}
-
-
-void  Pengrid::GetDofList(int* doflist,int* pnumberofdofspernode){
-
-	int j;
-	int doflist_per_node[MAXDOFSPERNODE];
-	int numberofdofspernode;
-	
-	node->GetDofList(&doflist_per_node[0],&numberofdofspernode);
-	for(j=0;j<numberofdofspernode;j++){
-		doflist[j]=doflist_per_node[j];
-	}
-
-	/*Assign output pointers:*/
-	*pnumberofdofspernode=numberofdofspernode;
-}
-
-void  Pengrid::PenaltyCreatePVectorThermal(Vec pg, void* vinputs, double kmax,int analysis_type,int sub_analysis_type){
-
-	const int numgrids=1;
-	const int NDOF1=1;
-	const int numdof=numgrids*NDOF1;
-	int       doflist[numdof];
-	double  P_terms[numdof]={0.0};
-	int    numberofdofspernode;
-	int    found=0;
-	double pressure;
-	int dofs1[1]={0};
-	double meltingpoint;
-	double beta;
-	double t_pmp;
-
-	ParameterInputs* inputs=NULL;
-
-	/*recover pointers: */
-	inputs=(ParameterInputs*)vinputs;
-
-	if(!active)return;
-
-	/*Get dof list: */
-	GetDofList(&doflist[0],&numberofdofspernode);
-
-	//First recover pressure
-	found=inputs->Recover("pressure",&pressure,1,dofs1,numgrids,(void**)&node);
-	if(!found)throw ErrorException(__FUNCT__," could not find pressure in inputs!");
-
-	//Compute pressure melting point
-	meltingpoint=matpar->GetMeltingPoint();
-	beta=matpar->GetBeta();
-	t_pmp=meltingpoint-beta*pressure;
-
-	//Add penalty load
-	P_terms[0]=kmax*pow((double)10,penalty_offset)*t_pmp;
-
-	/*Add P_terms to global vector pg: */
-	VecSetValues(pg,numdof,doflist,(const double*)P_terms,ADD_VALUES);
-}
-
+/*}}}1*/
+/*FUNCTION Pengrid::PenaltyCreatePVectorMelting {{{1*/
 void  Pengrid::PenaltyCreatePVectorMelting(Vec pg, void* vinputs, double kmax,int analysis_type,int sub_analysis_type){
 	
@@ -517,112 +612,52 @@
 	VecSetValues(pg,numdof,doflist,(const double*)P_terms,ADD_VALUES);
 }
-		
-
-#undef __FUNCT__ 
-#define __FUNCT__ "Pengrid::PenaltyConstrain"
-void  Pengrid::PenaltyConstrain(int* punstable,void* vinputs,int analysis_type,int sub_analysis_type){
-
-	if ((analysis_type==DiagnosticAnalysisEnum()) && ((sub_analysis_type==StokesAnalysisEnum()))){
-
-		/*No penalty to check*/
-		return;
-
-	}
-	else if (analysis_type==ThermalAnalysisEnum()){
-		
-		PenaltyConstrainThermal(punstable,vinputs,analysis_type,sub_analysis_type);
-		
-	}
-	else if (analysis_type==MeltingAnalysisEnum()){
-			
-		/*No penalty to check*/
-		return;
-
-	}
-	else{
-		throw ErrorException(__FUNCT__,exprintf("%s%i%s%i%s","analysis: ",analysis_type," and sub_analysis_type: ",sub_analysis_type," not supported yet"));
-	}
-
-}
-
-#undef __FUNCT__ 
-#define __FUNCT__ "Pengrid::PenaltyConstrainThermal"
-void  Pengrid::PenaltyConstrainThermal(int* punstable,void* vinputs,int analysis_type,int sub_analysis_type){
-
-	//   The penalty is stable if it doesn't change during to successive iterations.   
-
-	int found=0;
+/*}}}1*/
+/*FUNCTION Pengrid::PenaltyCreatePVectorThermal {{{1*/
+void  Pengrid::PenaltyCreatePVectorThermal(Vec pg, void* vinputs, double kmax,int analysis_type,int sub_analysis_type){
+
 	const int numgrids=1;
-
-
+	const int NDOF1=1;
+	const int numdof=numgrids*NDOF1;
+	int       doflist[numdof];
+	double  P_terms[numdof]={0.0};
+	int    numberofdofspernode;
+	int    found=0;
 	double pressure;
-	double temperature;
-	double beta,t_pmp;
+	int dofs1[1]={0};
 	double meltingpoint;
-	int    new_active;
-	int  dofs1[1]={0};
-	int  unstable=0;
-	int  reset_penalties=0;
-	
+	double beta;
+	double t_pmp;
+
 	ParameterInputs* inputs=NULL;
-
-	/*check that pengrid is not a clone (penalty to be added only once)*/
-	if (node->IsClone()){
-		unstable=0;
-		*punstable=unstable;
-		return;
-	}
 
 	/*recover pointers: */
 	inputs=(ParameterInputs*)vinputs;
 
-	//First recover beta, pressure and temperature vectors;
+	if(!active)return;
+
+	/*Get dof list: */
+	GetDofList(&doflist[0],&numberofdofspernode);
+
+	//First recover pressure
 	found=inputs->Recover("pressure",&pressure,1,dofs1,numgrids,(void**)&node);
 	if(!found)throw ErrorException(__FUNCT__," could not find pressure in inputs!");
-
-	found=inputs->Recover("temperature",&temperature,1,dofs1,numgrids,(void**)&node);
-	if(!found)throw ErrorException(__FUNCT__," could not find temperature in inputs!");
-	
-	found=inputs->Recover("reset_penalties",&reset_penalties);
-
-	if(reset_penalties)zigzag_counter=0;
 
 	//Compute pressure melting point
 	meltingpoint=matpar->GetMeltingPoint();
 	beta=matpar->GetBeta();
-
 	t_pmp=meltingpoint-beta*pressure;
 
-	//Figure out if temperature is over melting_point, in which case, this penalty needs to be activated.
-
-	if (temperature>t_pmp){
-		new_active=1;
-	}
-	else{
-		new_active=0;
-	}
-
-
-	//Figure out stability of this penalty
-	if (active==new_active){
-		unstable=0;
-	}
-	else{
-		unstable=1;
-		if(stabilize_constraints)zigzag_counter++;
-	}
-
-	/*If penalty keeps zigzagging more than 5 times: */
-	if(stabilize_constraints){
-		if(zigzag_counter>stabilize_constraints){
-			unstable=0;
-			active=1;
-		}
-	}
-
-	//Set penalty flag
-	active=new_active;
-
-	//*Assign output pointers:*/
-	*punstable=unstable;
-}
+	//Add penalty load
+	P_terms[0]=kmax*pow((double)10,penalty_offset)*t_pmp;
+
+	/*Add P_terms to global vector pg: */
+	VecSetValues(pg,numdof,doflist,(const double*)P_terms,ADD_VALUES);
+}
+/*}}}1*/
+/*FUNCTION Pengrid::UpdateFromInputs {{{1*/
+#undef __FUNCT__ 
+#define __FUNCT__ "Pengrid::UpdateFromInputs"
+void  Pengrid::UpdateFromInputs(void* inputs){
+	
+}
+/*}}}1*/
Index: /issm/trunk/src/c/objects/Penpair.cpp
===================================================================
--- /issm/trunk/src/c/objects/Penpair.cpp	(revision 2906)
+++ /issm/trunk/src/c/objects/Penpair.cpp	(revision 2907)
@@ -18,8 +18,11 @@
 
 		
+/*Object constructors and destructor*/
+/*FUNCTION Penpair::constructor {{{1*/
 Penpair::Penpair(){
 	return;
 }
-
+/*}}}1*/
+/*FUNCTION Penpair::creation {{{1*/
 Penpair::Penpair(int	penpair_id, double penpair_penalty_offset,int  penpair_penalty_lock,int penpair_node_ids[2],int penpair_dof){
 	
@@ -39,10 +42,126 @@
 	return;
 }
-
+/*}}}1*/
+/*FUNCTION Penpair::destructor {{{1*/
 Penpair::~Penpair(){
 	return;
 }
+/*}}}1*/
 		
-void Penpair::Echo(void){
+/*Object marshall*/
+/*FUNCTION Penpair::Marshall {{{1*/
+void  Penpair::Marshall(char** pmarshalled_dataset){
+
+	char* marshalled_dataset=NULL;
+	int   enum_type=0;
+
+	/*recover marshalled_dataset: */
+	marshalled_dataset=*pmarshalled_dataset;
+
+	/*get enum type of Penpair: */
+	enum_type=PenpairEnum();
+	
+	/*marshall enum: */
+	memcpy(marshalled_dataset,&enum_type,sizeof(enum_type));marshalled_dataset+=sizeof(enum_type);
+	
+	/*marshall Penpair data: */
+	memcpy(marshalled_dataset,&id,sizeof(id));marshalled_dataset+=sizeof(id);
+	memcpy(marshalled_dataset,&penalty_offset,sizeof(penalty_offset));marshalled_dataset+=sizeof(penalty_offset);
+	memcpy(marshalled_dataset,&penalty_lock,sizeof(penalty_lock));marshalled_dataset+=sizeof(penalty_lock);
+	memcpy(marshalled_dataset,&dof,sizeof(dof));marshalled_dataset+=sizeof(dof);
+	memcpy(marshalled_dataset,&node_ids,sizeof(node_ids));marshalled_dataset+=sizeof(node_ids);
+	memcpy(marshalled_dataset,&node_offsets,sizeof(node_offsets));marshalled_dataset+=sizeof(node_offsets);
+
+	*pmarshalled_dataset=marshalled_dataset;
+	return;
+}
+/*}}}1*/
+/*FUNCTION Penpair::MarshallSize {{{1*/
+int   Penpair::MarshallSize(){
+
+	return sizeof(id)+
+		sizeof(penalty_offset)+
+		sizeof(penalty_lock)+
+		sizeof(dof)+
+		sizeof(node_ids)+
+		sizeof(node_offsets)+
+		sizeof(int); //sizeof(int) for enum type
+}
+/*}}}1*/
+/*FUNCTION Penpair::Demarshall {{{1*/
+void  Penpair::Demarshall(char** pmarshalled_dataset){
+
+	int i;
+	char* marshalled_dataset=NULL;
+
+	/*recover marshalled_dataset: */
+	marshalled_dataset=*pmarshalled_dataset;
+
+	/*this time, no need to get enum type, the pointer directly points to the beginning of the 
+	 *object data (thanks to DataSet::Demarshall):*/
+
+	memcpy(&id,marshalled_dataset,sizeof(id));marshalled_dataset+=sizeof(id);
+	memcpy(&penalty_offset,marshalled_dataset,sizeof(penalty_offset));marshalled_dataset+=sizeof(penalty_offset);
+	memcpy(&penalty_lock,marshalled_dataset,sizeof(penalty_lock));marshalled_dataset+=sizeof(penalty_lock);
+	memcpy(&dof,marshalled_dataset,sizeof(dof));marshalled_dataset+=sizeof(dof);
+	memcpy(&node_ids,marshalled_dataset,sizeof(node_ids));marshalled_dataset+=sizeof(node_ids);
+	memcpy(&node_offsets,marshalled_dataset,sizeof(node_offsets));marshalled_dataset+=sizeof(node_offsets);
+
+	for(i=0;i<2;i++){
+		nodes[i]=NULL;
+	}
+
+	/*return: */
+	*pmarshalled_dataset=marshalled_dataset;
+	return;
+}
+/*}}}1*/
+
+/*Object functions*/
+/*FUNCTION Penpair::Configure {{{1*/
+#undef __FUNCT__
+#define __FUNCT__ "Penpair::Configure"
+
+void  Penpair::Configure(void* pelementsin,void* pnodesin,void* pmaterialsin){
+
+	DataSet* elementsin=NULL;
+	DataSet* nodesin=NULL;
+
+	/*Recover pointers :*/
+	elementsin=(DataSet*)pelementsin;
+	nodesin=(DataSet*)pnodesin;
+	/*Link this load with its nodes: */
+	ResolvePointers((Object**)nodes,node_ids,node_offsets,2,nodesin);
+
+}
+/*}}}1*/
+/*FUNCTION Penpair::copy {{{1*/
+Object* Penpair::copy() {
+	return new Penpair(*this); 
+}
+/*}}}1*/
+/*FUNCTION Penpair::CreateKMatrix {{{1*/
+#undef __FUNCT__ 
+#define __FUNCT__ "Penpair::CreateKMatrix"
+
+void  Penpair::CreateKMatrix(Mat Kgg,void* inputs,int analysis_type,int sub_analysis_type){
+
+	/*No loads applied, do nothing: */
+	return;
+
+}
+/*}}}1*/
+/*FUNCTION Penpair::CreatePVector {{{1*/
+#undef __FUNCT__ 
+#define __FUNCT__ "Penpair::CreatePVector"
+void  Penpair::CreatePVector(Vec pg, void* inputs, int analysis_type,int sub_analysis_type){
+
+	/*No loads applied, do nothing: */
+	return;
+
+}
+/*}}}1*/
+/*FUNCTION Penpair::DeepEcho {{{1*/
+void Penpair::DeepEcho(void){
 
 	int i;
@@ -60,6 +179,8 @@
 	}
 	return;
-}
-void Penpair::DeepEcho(void){
+}		
+/*}}}1*/
+/*FUNCTION Penpair::Echo {{{1*/
+void Penpair::Echo(void){
 
 	int i;
@@ -77,125 +198,27 @@
 	}
 	return;
-}		
-void  Penpair::Marshall(char** pmarshalled_dataset){
-
-	char* marshalled_dataset=NULL;
-	int   enum_type=0;
-
-	/*recover marshalled_dataset: */
-	marshalled_dataset=*pmarshalled_dataset;
-
-	/*get enum type of Penpair: */
-	enum_type=PenpairEnum();
-	
-	/*marshall enum: */
-	memcpy(marshalled_dataset,&enum_type,sizeof(enum_type));marshalled_dataset+=sizeof(enum_type);
-	
-	/*marshall Penpair data: */
-	memcpy(marshalled_dataset,&id,sizeof(id));marshalled_dataset+=sizeof(id);
-	memcpy(marshalled_dataset,&penalty_offset,sizeof(penalty_offset));marshalled_dataset+=sizeof(penalty_offset);
-	memcpy(marshalled_dataset,&penalty_lock,sizeof(penalty_lock));marshalled_dataset+=sizeof(penalty_lock);
-	memcpy(marshalled_dataset,&dof,sizeof(dof));marshalled_dataset+=sizeof(dof);
-	memcpy(marshalled_dataset,&node_ids,sizeof(node_ids));marshalled_dataset+=sizeof(node_ids);
-	memcpy(marshalled_dataset,&node_offsets,sizeof(node_offsets));marshalled_dataset+=sizeof(node_offsets);
-
-	*pmarshalled_dataset=marshalled_dataset;
-	return;
-}
-		
-int   Penpair::MarshallSize(){
-
-	return sizeof(id)+
-		sizeof(penalty_offset)+
-		sizeof(penalty_lock)+
-		sizeof(dof)+
-		sizeof(node_ids)+
-		sizeof(node_offsets)+
-		sizeof(int); //sizeof(int) for enum type
-}
-
+}
+/*}}}1*/
+/*FUNCTION Penpair::Enum {{{1*/
+int Penpair::Enum(void){
+
+	return PenpairEnum();
+}
+/*}}}1*/
+/*FUNCTION Penpair::GetId {{{1*/
+int    Penpair::GetId(void){ return id; }
+/*}}}1*/
+/*FUNCTION Penpair::GetName {{{1*/
 char* Penpair::GetName(void){
 	return "penpair";
 }
-		
-
-void  Penpair::Demarshall(char** pmarshalled_dataset){
-
-	int i;
-	char* marshalled_dataset=NULL;
-
-	/*recover marshalled_dataset: */
-	marshalled_dataset=*pmarshalled_dataset;
-
-	/*this time, no need to get enum type, the pointer directly points to the beginning of the 
-	 *object data (thanks to DataSet::Demarshall):*/
-
-	memcpy(&id,marshalled_dataset,sizeof(id));marshalled_dataset+=sizeof(id);
-	memcpy(&penalty_offset,marshalled_dataset,sizeof(penalty_offset));marshalled_dataset+=sizeof(penalty_offset);
-	memcpy(&penalty_lock,marshalled_dataset,sizeof(penalty_lock));marshalled_dataset+=sizeof(penalty_lock);
-	memcpy(&dof,marshalled_dataset,sizeof(dof));marshalled_dataset+=sizeof(dof);
-	memcpy(&node_ids,marshalled_dataset,sizeof(node_ids));marshalled_dataset+=sizeof(node_ids);
-	memcpy(&node_offsets,marshalled_dataset,sizeof(node_offsets));marshalled_dataset+=sizeof(node_offsets);
-
-	for(i=0;i<2;i++){
-		nodes[i]=NULL;
-	}
-
-	/*return: */
-	*pmarshalled_dataset=marshalled_dataset;
-	return;
-}
-int Penpair::Enum(void){
-
-	return PenpairEnum();
-}
-
-int    Penpair::GetId(void){ return id; }
-
+/*}}}1*/
+/*FUNCTION Penpair::MyRank {{{1*/
 int    Penpair::MyRank(void){ 
 	extern int my_rank;
 	return my_rank; 
 }
-
-#undef __FUNCT__
-#define __FUNCT__ "Penpair::Configure"
-
-void  Penpair::Configure(void* pelementsin,void* pnodesin,void* pmaterialsin){
-
-	DataSet* elementsin=NULL;
-	DataSet* nodesin=NULL;
-
-	/*Recover pointers :*/
-	elementsin=(DataSet*)pelementsin;
-	nodesin=(DataSet*)pnodesin;
-	/*Link this load with its nodes: */
-	ResolvePointers((Object**)nodes,node_ids,node_offsets,2,nodesin);
-
-}
-
-#undef __FUNCT__ 
-#define __FUNCT__ "Penpair::CreateKMatrix"
-
-void  Penpair::CreateKMatrix(Mat Kgg,void* inputs,int analysis_type,int sub_analysis_type){
-
-	/*No loads applied, do nothing: */
-	return;
-
-}
-
-#undef __FUNCT__ 
-#define __FUNCT__ "Penpair::CreatePVector"
-void  Penpair::CreatePVector(Vec pg, void* inputs, int analysis_type,int sub_analysis_type){
-
-	/*No loads applied, do nothing: */
-	return;
-
-}
-#undef __FUNCT__ 
-#define __FUNCT__ "Penpair::UpdateFromInputs"
-void  Penpair::UpdateFromInputs(void* inputs){
-	
-}
-
+/*}}}1*/
+/*FUNCTION Penpair::PenaltyCreateKMatrix {{{1*/
 #undef __FUNCT__ 
 #define __FUNCT__ "Penpair::PenaltyCreateKMatrix"
@@ -205,5 +228,6 @@
 	return;
 }
-		
+/*}}}1*/
+/*FUNCTION Penpair::PenaltyCreatePVector {{{1*/
 #undef __FUNCT__ 
 #define __FUNCT__ "Penpair::PenaltyCreatePVector"
@@ -212,7 +236,10 @@
 	return;
 }
-
-Object* Penpair::copy() {
-	return new Penpair(*this); 
-}
-
+/*}}}1*/
+/*FUNCTION Penpair::UpdateFromInputs {{{1*/
+#undef __FUNCT__ 
+#define __FUNCT__ "Penpair::UpdateFromInputs"
+void  Penpair::UpdateFromInputs(void* inputs){
+	
+}
+/*}}}1*/
Index: /issm/trunk/src/c/objects/Result.cpp
===================================================================
--- /issm/trunk/src/c/objects/Result.cpp	(revision 2906)
+++ /issm/trunk/src/c/objects/Result.cpp	(revision 2907)
@@ -18,8 +18,19 @@
 #include "../include/typedefs.h"
 
+/*Object constructors and destructor*/
+/*FUNCTION Result::constructor {{{1*/
 Result::Result(){
 	return;
 }
-
+/*}}}1*/
+/*FUNCTION Result::destructor {{{1*/
+Result::~Result(){
+	xfree((void**)&fieldname);
+	VecFree(&field);
+	xfree((void**)&dfield);
+	xfree((void**)&cfield);
+}
+/*}}}1*/
+/*FUNCTION Result::Result(const Result& result) {{{1*/
 Result::Result(const Result& result){
 
@@ -52,5 +63,6 @@
 	}
 }
-
+/*}}}1*/
+/*FUNCTION Result::Result(int result_id,double result_time,int result_step,char* result_fieldname,Vec result_field) {{{1*/
 #undef __FUNCT__ 
 #define __FUNCT__ "Result::Result"
@@ -73,5 +85,6 @@
 	cfield=NULL;
 }
-
+/*}}}1*/
+/*FUNCTION Result::Result(int result_id,double result_time,int result_step,char* result_fldname,double* result_field,int result_size){{{1*/
 #undef __FUNCT__ 
 #define __FUNCT__ "Result::Result"
@@ -96,5 +109,6 @@
 	cfield=NULL;
 }
-
+/*}}}1*/
+/*FUNCTION Result::Result(int result_id,double result_time,int result_step,char* result_fieldname,char* result_field) {{{1*/
 #undef __FUNCT__ 
 #define __FUNCT__ "Result::Result"
@@ -118,13 +132,41 @@
 	dfield=NULL;
 }
-
-Result::~Result(){
-	xfree((void**)&fieldname);
-	VecFree(&field);
-	xfree((void**)&dfield);
-	xfree((void**)&cfield);
-}
+/*}}}1*/
 		
-void Result::Echo(void){
+/*Object marshall*/
+/*FUNCTION Result::Marshall {{{1*/
+#undef __FUNCT__ 
+#define __FUNCT__ "Result::Marshall"
+void  Result::Marshall(char** pmarshalled_dataset){
+
+	throw  ErrorException(__FUNCT__," not supported yet!");
+
+}
+/*}}}1*/
+/*FUNCTION Result::MarshallSize {{{1*/
+#undef __FUNCT__ 
+#define __FUNCT__ "Result::MarshallSize"
+int   Result::MarshallSize(){
+
+	throw  ErrorException(__FUNCT__," not supported yet!");
+}
+/*}}}1*/
+/*FUNCTION Result::Demarshall {{{1*/
+#undef __FUNCT__ 
+#define __FUNCT__ "Result::Demarshall"
+void  Result::Demarshall(char** pmarshalled_dataset){
+
+	throw  ErrorException(__FUNCT__," not supported yet!");
+}
+/*}}}1*/
+
+/*Object functions*/
+/*FUNCTION Result::copy {{{1*/
+Object* Result::copy() {
+	return new Result(*this); 
+}
+/*}}}1*/
+/*FUNCTION Result::DeepEcho {{{1*/
+void Result::DeepEcho(void){
 
 	printf("Result:\n");
@@ -144,7 +186,8 @@
 		printf("   field string %s\n",cfield);
 	}
-}
-
-void Result::DeepEcho(void){
+}	
+/*}}}1*/
+/*FUNCTION Result::Echo {{{1*/
+void Result::Echo(void){
 
 	printf("Result:\n");
@@ -164,40 +207,61 @@
 		printf("   field string %s\n",cfield);
 	}
-}	
-
-#undef __FUNCT__ 
-#define __FUNCT__ "Result::Marshall"
-void  Result::Marshall(char** pmarshalled_dataset){
-
-	throw  ErrorException(__FUNCT__," not supported yet!");
-
-}
-	
-#undef __FUNCT__ 
-#define __FUNCT__ "Result::MarshallSize"
-int   Result::MarshallSize(){
-
-	throw  ErrorException(__FUNCT__," not supported yet!");
-}
-
+}
+/*}}}1*/
+/*FUNCTION Result::Enum {{{1*/
+int Result::Enum(void){
+
+	return ResultEnum();
+
+}
+/*}}}1*/
+/*FUNCTION Result::GetFieldName {{{1*/
+char*    Result::GetFieldName(){
+	return fieldname;
+}
+/*}}}1*/
+/*FUNCTION Result::GetField(Vec* pfield) {{{1*/
+void  Result::GetField(Vec* pfield){
+
+	VecDuplicatePatch(pfield,field);
+
+}
+/*}}}1*/
+/*FUNCTION Result::GetField(char** pcfield) {{{1*/
+void  Result::GetField(char** pcfield){
+
+	char* string=NULL;
+
+	string=(char*)xmalloc((strlen(cfield)+1)*sizeof(char));
+	strcpy(string,cfield);
+
+	*pcfield=string;
+}
+/*}}}1*/
+/*FUNCTION Result::GetField(double** pfield) {{{1*/
+void  Result::GetField(double** pfield){
+	*pfield=(double*)xmalloc(size*sizeof(double));
+	memcpy(*pfield,dfield,size*sizeof(double));
+}
+/*}}}1*/
+/*FUNCTION Result::GetId {{{1*/
+int    Result::GetId(void){ return id; }
+/*}}}1*/
+/*FUNCTION Result::GetName {{{1*/
 char* Result::GetName(void){
 	return "result";
 }
-		
-
-#undef __FUNCT__ 
-#define __FUNCT__ "Result::Demarshall"
-void  Result::Demarshall(char** pmarshalled_dataset){
-
-	throw  ErrorException(__FUNCT__," not supported yet!");
-}
-
-int Result::Enum(void){
-
-	return ResultEnum();
-
-}
-int    Result::GetId(void){ return id; }
-
+/*}}}1*/
+/*FUNCTION Result::GetStep {{{1*/
+int    Result::GetStep(){
+	return step;
+}
+/*}}}1*/
+/*FUNCTION Result::GetTime {{{1*/
+double Result::GetTime(){
+	return time;
+}
+/*}}}1*/
+/*FUNCTION Result::MyRank {{{1*/
 int    Result::MyRank(void){ 
 	extern int my_rank;
@@ -205,17 +269,6 @@
 	return my_rank; 
 }
-
-double Result::GetTime(){
-	return time;
-}
-
-int    Result::GetStep(){
-	return step;
-}
-
-char*    Result::GetFieldName(){
-	return fieldname;
-}
-		
+/*}}}1*/
+/*FUNCTION Result::WriteData {{{1*/
 void   Result::WriteData(FILE* fid){
 
@@ -247,30 +300,3 @@
 	}
 }
-
-
-void  Result::GetField(Vec* pfield){
-
-	VecDuplicatePatch(pfield,field);
-
-}
-
-void  Result::GetField(char** pcfield){
-
-	char* string=NULL;
-
-	string=(char*)xmalloc((strlen(cfield)+1)*sizeof(char));
-	strcpy(string,cfield);
-
-	*pcfield=string;
-}
-
-void  Result::GetField(double** pfield){
-	*pfield=(double*)xmalloc(size*sizeof(double));
-	memcpy(*pfield,dfield,size*sizeof(double));
-}
-
-Object* Result::copy() {
-	return new Result(*this); 
-}
-		
-
+/*}}}1*/
Index: /issm/trunk/src/c/objects/Rgb.cpp
===================================================================
--- /issm/trunk/src/c/objects/Rgb.cpp	(revision 2906)
+++ /issm/trunk/src/c/objects/Rgb.cpp	(revision 2907)
@@ -16,8 +16,11 @@
 #include "./objects.h"
 
-		
+/*Object constructors and destructor*/
+/*FUNCTION Rgb::constructor {{{1*/
 Rgb::Rgb(){
 	return;
 }
+/*}}}1*/
+/*FUNCTION Rgb::creation {{{1*/
 Rgb::Rgb(int rgb_id,int rgb_nodeid1,int rgb_nodeid2, int rgb_dof){
 
@@ -29,27 +32,13 @@
 	return;
 }
-
+/*}}}1*/
+/*FUNCTION Rgb::destructor {{{1*/
 Rgb::~Rgb(){
 	return;
 }
+/*}}}1*/
 		
-void Rgb::Echo(void){
-
-	printf("Rgb:\n");
-	printf("   id: %i\n",id);
-	printf("   nodeid1: %i\n",nodeid1);
-	printf("   nodeid2: %i\n",nodeid2);
-	printf("   dof: %i\n",dof);
-	return;
-}
-void Rgb::DeepEcho(void){
-
-	printf("Rgb:\n");
-	printf("   id: %i\n",id);
-	printf("   nodeid1: %i\n",nodeid1);
-	printf("   nodeid2: %i\n",nodeid2);
-	printf("   dof: %i\n",dof);
-	return;
-}		
+/*Object marshall*/
+/*FUNCTION Rgb::Marshall {{{1*/
 void  Rgb::Marshall(char** pmarshalled_dataset){
 
@@ -75,5 +64,6 @@
 	return;
 }
-		
+/*}}}1*/
+/*FUNCTION Rgb::MarshallSize {{{1*/
 int   Rgb::MarshallSize(){
 
@@ -84,10 +74,6 @@
 		sizeof(int); //sizeof(int) for enum type
 }
-
-char* Rgb::GetName(void){
-	return "rgb";
-}
-		
-
+/*}}}1*/
+/*FUNCTION Rgb::Demarshall {{{1*/
 void  Rgb::Demarshall(char** pmarshalled_dataset){
 
@@ -109,4 +95,35 @@
 	return;
 }
+/*}}}1*/
+
+/*Object functions*/
+/*FUNCTION Rgb::copy {{{1*/
+Object* Rgb::copy() {
+	return new Rgb(*this); 
+}
+/*}}}1*/
+/*FUNCTION Rgb::DeepEcho {{{1*/
+void Rgb::DeepEcho(void){
+
+	printf("Rgb:\n");
+	printf("   id: %i\n",id);
+	printf("   nodeid1: %i\n",nodeid1);
+	printf("   nodeid2: %i\n",nodeid2);
+	printf("   dof: %i\n",dof);
+	return;
+}		
+/*}}}1*/
+/*FUNCTION Rgb::Echo {{{1*/
+void Rgb::Echo(void){
+
+	printf("Rgb:\n");
+	printf("   id: %i\n",id);
+	printf("   nodeid1: %i\n",nodeid1);
+	printf("   nodeid2: %i\n",nodeid2);
+	printf("   dof: %i\n",dof);
+	return;
+}
+/*}}}1*/
+/*FUNCTION Rgb::Enum {{{1*/
 int Rgb::Enum(void){
 
@@ -114,25 +131,34 @@
 
 }
+/*}}}1*/
+/*FUNCTION Rgb::GetId {{{1*/
 int    Rgb::GetId(void){ return id; }
-
+/*}}}1*/
+/*FUNCTION Rgb::GetName {{{1*/
+char* Rgb::GetName(void){
+	return "rgb";
+}
+/*}}}1*/
+/*FUNCTION Rgb::GetNodeId1{{{1*/
+int   Rgb::GetNodeId1(){
+	
+	return nodeid1;
+}
+/*}}}1*/
+/*FUNCTION Rgb::GetNodeId2 {{{1*/
+int   Rgb::GetNodeId2(){
+	
+	return nodeid2;
+}
+/*}}}1*/
+/*FUNCTION Rgb::GetDof {{{1*/
+int Rgb::GetDof(){
+	return dof;
+}
+/*}}}1*/
+/*FUNCTION Rgb::MyRank {{{1*/
 int    Rgb::MyRank(void){ 
 	extern int my_rank;
 	return my_rank; 
 }
-
-int   Rgb::GetNodeId1(){
-	
-	return nodeid1;
-}
-int   Rgb::GetNodeId2(){
-	
-	return nodeid2;
-}
-
-int Rgb::GetDof(){
-	return dof;
-}
-Object* Rgb::copy() {
-	return new Rgb(*this); 
-}
-
+/*}}}1*/