Index: issm/trunk/src/c/objects/Node.cpp
===================================================================
--- issm/trunk/src/c/objects/Node.cpp	(revision 832)
+++ issm/trunk/src/c/objects/Node.cpp	(revision 848)
@@ -550,6 +550,6 @@
 		
 #undef __FUNCT__ 
-#define __FUNCT__ "Node::VelocityDepthAverageAtBase"
-void  Node::VelocityDepthAverageAtBase(Vec ug,double* ug_serial){
+#define __FUNCT__ "Node::FieldDepthAverageAtBase"
+void  Node::FieldDepthAverageAtBase(Vec field,double* field_serial,char* fieldname){
 
 	/* node data: */
@@ -560,11 +560,8 @@
 	Node* node=NULL;
 	Node* upper_node=NULL;
-	double velocity2[2];
-	double velocity1[2];
-	double velocity_average[2];
-	double sum[2];
 	double z1,z2,dz;
 	double thickness;
 
+	/*Are we on the base, not on the surface, and not on a clone node?:*/
 	
 	if(onbed==1 & clone==0 &onsurface==0){
@@ -572,53 +569,206 @@
 		doflist1=GetDofList1();
 
-		/*this node is on the bed. We are going to:
-		 * follow the upper nodes until we reach the surface. At each upper node, we'll grab the 
-		 * velocity for this node, and add it to ug. We'll finish by 
-		 * we grab the velocity. Once we know the velocity, we follow the upper nodes, 
-		 * inserting the same velocity value into ug, until we reach the surface: */
-		sum[0]=0;
-		sum[1]=0;
-		thickness=0;
-
-		/*get dofs for this base node velocity: we know there are two dofs in ug_serial */
-		dofx=2*doflist1;
-		dofy=2*doflist1+1;
-
-		node=this;
-		for(;;){
-
-			if (node->IsOnSurface())break;
-
-			doflist1=node->GetDofList1();
+		/*this node is on the bed. We are going to, follow the upper nodes until we reach the surface. At each upper node, 
+		 * we'll grab the * field for this node, and add it to overall field: */
+
+		if(strcmp(fieldname,"velocity")==0){
+
+			/*field is a velocity, 2 dofs per node: */
+			double velocity2[2];
+			double velocity1[2];
+			double velocity_average[2];
+			double sum[2];
+
+			sum[0]=0;
+			sum[1]=0;
+			thickness=0;
+
+			/*get dofs for this base node velocity: we know there are two dofs in field_serial */
+			dofx=2*doflist1;
+			dofy=2*doflist1+1;
+
+			node=this;
+			for(;;){
+
+				if (node->IsOnSurface())break;
+
+				doflist1=node->GetDofList1();
+				
+				velocity1[0]=field_serial[2*doflist1];
+				velocity1[1]=field_serial[2*doflist1+1];
+				z1=node->GetZ();
+
+				upper_node=node->GetUpperNode();
+				doflist1=upper_node->GetDofList1();
 			
-			velocity1[0]=ug_serial[2*doflist1];
-			velocity1[1]=ug_serial[2*doflist1+1];
-			z1=node->GetZ();
-
-			upper_node=node->GetUpperNode();
-			doflist1=upper_node->GetDofList1();
-		
-			velocity2[0]=ug_serial[2*doflist1];
-			velocity2[1]=ug_serial[2*doflist1+1];
-			z2=upper_node->GetZ();
-
-			dz=(z2-z1);
-			thickness+=dz;
-			velocity_average[0]=(velocity1[0]+velocity2[0])/2.0;
-			velocity_average[1]=(velocity1[1]+velocity2[1])/2.0;
-
-			sum[0]+=velocity_average[0]*dz;
-			sum[1]+=velocity_average[1]*dz;
-
-			/* get next node: */
-			node=node->GetUpperNode();
+				velocity2[0]=field_serial[2*doflist1];
+				velocity2[1]=field_serial[2*doflist1+1];
+				z2=upper_node->GetZ();
+
+				dz=(z2-z1);
+				thickness+=dz;
+				velocity_average[0]=(velocity1[0]+velocity2[0])/2.0;
+				velocity_average[1]=(velocity1[1]+velocity2[1])/2.0;
+
+				sum[0]+=velocity_average[0]*dz;
+				sum[1]+=velocity_average[1]*dz;
+
+				/* get next node: */
+				node=node->GetUpperNode();
+			}
+
+			sum[0]=sum[0]/thickness;
+			sum[1]=sum[1]/thickness;
+
+			/* Plfield velocity_average*deltaH/H into base of field: */
+			VecSetValues(field,1,&dofx,&sum[0],INSERT_VALUES);
+			VecSetValues(field,1,&dofy,&sum[1],INSERT_VALUES);
 		}
-
-		sum[0]=sum[0]/thickness;
-		sum[1]=sum[1]/thickness;
-
-		/* Plug velocity_average*deltaH/H into base of ug: */
-		VecSetValues(ug,1,&dofx,&sum[0],INSERT_VALUES);
-		VecSetValues(ug,1,&dofy,&sum[1],INSERT_VALUES);
-	}
-}
+		else{
+			/*field is regular, 1 dof per node: */
+			double field2;
+			double field1;
+			double field_average;
+			double sum;
+
+			sum=0;
+			thickness=0;
+
+			/*get dofs for this base node velocity: we know there are two dofs in field_serial */
+			dofx=doflist1;
+
+			node=this;
+			for(;;){
+
+				if (node->IsOnSurface())break;
+
+				doflist1=node->GetDofList1();
+				
+				field1=field_serial[doflist1];
+				z1=node->GetZ();
+
+				upper_node=node->GetUpperNode();
+				doflist1=upper_node->GetDofList1();
+			
+				field2=field_serial[2*doflist1];
+				z2=upper_node->GetZ();
+
+				dz=(z2-z1);
+				thickness+=dz;
+				field_average=(field1+field2)/2.0;
+
+				sum+=field_average*dz;
+
+				/* get next node: */
+				node=node->GetUpperNode();
+			}
+
+			sum=sum/thickness;
+
+			/* Plug field_average*deltH/H into base of field: */
+			VecSetValues(field,1,&dofx,&sum,INSERT_VALUES);
+		}
+	}
+}
+
+#undef __FUNCT__ 
+#define __FUNCT__ "Node::FieldExtrude"
+void  Node::FieldExtrude(Vec field,double* field_serial,char* field_name){
+		
+	/* node data: */
+	int   numberofdofspernode;
+	Node* node=NULL;
+	int   i;
+
+	/*Is this node on bed? :*/
+	if (onbed){
+
+		if (strcmp(field_name,"velocity")==0){
+
+			/* node data: */
+			const int    numdof=2;
+			int          doflist[numdof];
+			int          nodedofs[2];
+			double       fieldnode[2];
+
+
+			GetDofList(&doflist[0],&numberofdofspernode);
+
+			//initilaize node
+			node=this;
+		
+			/*get field for this base node: */
+			fieldnode[0]=field_serial[doflist[0]];
+			fieldnode[1]=field_serial[doflist[1]];
+
+			//go through all nodes which sit on top of this node, until we reach the surface, 
+			//and plfield  field in field
+			for(;;){
+
+				node->GetDofList(&nodedofs[0],&numberofdofspernode);
+				VecSetValues(field,1,&nodedofs[0],&fieldnode[0],INSERT_VALUES);
+				VecSetValues(field,1,&nodedofs[1],&fieldnode[1],INSERT_VALUES);
+
+				if (node->IsOnSurface())break;
+				/*get next node: */
+				node=node->GetUpperNode();
+			}
+		} //if (strcmp(field_name,"velocity")==0)
+		else if ( 
+				(strcmp(field_name,"thickness")==0) ||
+				(strcmp(field_name,"surface")==0)  ||
+				(strcmp(field_name,"bed")==0)  || 
+				(strcmp(field_name,"slopex")==0)  ||
+				(strcmp(field_name,"slopey")==0) 
+				){
+
+			/* node data: */
+			int          doflist;
+			int          nodedofs;
+			double       fieldnode;
+				
+			GetDofList(&doflist,&numberofdofspernode);
+
+			//initilaize node
+			node=this;
+
+			/*get field for this bed node: */
+			fieldnode=field_serial[doflist];
+
+			//go through all nodes which sit on top of this node, until we reach the surface, 
+			//and pltg  fieldnode in field:
+			for(;;){
+
+				node->GetDofList(&nodedofs,&numberofdofspernode);
+				VecSetValues(field,1,&nodedofs,&fieldnode,INSERT_VALUES);
+
+				if (node->IsOnSurface())break;
+				/*get next node: */
+				node=node->GetUpperNode();
+			}
+		}
+		else throw ErrorException(__FUNCT__,exprintf("%s%s%s"," field ",field_name," not supported yet!"));
+
+	} //if (extrude)
+}
+		
+void  Node::UpdateNodePosition(double* thickness,double* bed){
+
+	int dof;
+
+	dof=this->GetDofList1();
+
+	if(onbed){
+
+		/*this node is on the bed, set its z position to the new bed: */
+		this->x[2]=bed[dof];
+	}
+	else if (onsurface){
+
+		/*this node is on the surface, set its z position to the new bed+ new thickness: */
+		this->x[2]=bed[dof]+thickness[dof];
+	}
+	else{
+		/*do nothing, we only update grids belonging to the exterior layers: */
+	}
+}
