Index: /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp
===================================================================
--- /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp	(revision 25553)
+++ /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp	(revision 25554)
@@ -575,6 +575,6 @@
 			xDelete<int>(approximations);
 
-			for(int i=0;i<nodes->Size();i++){
-				Node* node=xDynamicCast<Node*>(nodes->GetObjectByOffset(i));
+			for(Object* & object: nodes->objects){
+				Node* node=xDynamicCast<Node*>(object);
 				int   sid = node->Sid();
 				if(sid>=iomodel->numberofvertices+iomodel->numberofelements){
@@ -975,10 +975,10 @@
 		case 3:
 			parameters->AddObject(iomodel->CopyConstantObject("md.friction.coupling",FrictionCouplingEnum));
-			parameters->AddObject(iomodel->CopyConstantObject("md.friction.effective_pressure_limit",FrictionEffectivePressureLimitEnum));			
+			parameters->AddObject(iomodel->CopyConstantObject("md.friction.effective_pressure_limit",FrictionEffectivePressureLimitEnum));
 			break;
 		case 4:
 			parameters->AddObject(iomodel->CopyConstantObject("md.friction.gamma",FrictionGammaEnum));
 			parameters->AddObject(iomodel->CopyConstantObject("md.friction.coupling",FrictionCouplingEnum));
-			parameters->AddObject(iomodel->CopyConstantObject("md.friction.effective_pressure_limit",FrictionEffectivePressureLimitEnum));			
+			parameters->AddObject(iomodel->CopyConstantObject("md.friction.effective_pressure_limit",FrictionEffectivePressureLimitEnum));
 			break;
 		case 5:
@@ -996,5 +996,5 @@
 		case 9:
 			parameters->AddObject(iomodel->CopyConstantObject("md.friction.gamma",FrictionGammaEnum));
-			parameters->AddObject(iomodel->CopyConstantObject("md.friction.effective_pressure_limit",FrictionEffectivePressureLimitEnum));			
+			parameters->AddObject(iomodel->CopyConstantObject("md.friction.effective_pressure_limit",FrictionEffectivePressureLimitEnum));
 			break;
 		case 10:
@@ -1008,9 +1008,9 @@
 		case 11:
 			parameters->AddObject(new IntParam(FrictionCouplingEnum,2));
-			parameters->AddObject(iomodel->CopyConstantObject("md.friction.effective_pressure_limit",FrictionEffectivePressureLimitEnum));			
+			parameters->AddObject(iomodel->CopyConstantObject("md.friction.effective_pressure_limit",FrictionEffectivePressureLimitEnum));
 			break;
 		case 12:
 			parameters->AddObject(new IntParam(FrictionCouplingEnum,2));
-			parameters->AddObject(iomodel->CopyConstantObject("md.friction.effective_pressure_limit",FrictionEffectivePressureLimitEnum));			
+			parameters->AddObject(iomodel->CopyConstantObject("md.friction.effective_pressure_limit",FrictionEffectivePressureLimitEnum));
 			break;
 		default: _error_("Friction law "<<frictionlaw<<" not implemented yet");
@@ -1049,5 +1049,5 @@
 		#endif
 
-		
+
 		if(is_schur_cg_solver)
 		 solutionsequence_schurcg(femmodel);
@@ -1187,5 +1187,5 @@
 			return;
 		default:
-			_error_("Approximation "<<EnumToStringx(approximation)<<" not supported");
+			_error_("Approximation "<<EnumToStringx(approximation)<<"("<<approximation<<") not supported");
 	}
 }/*}}}*/
@@ -1721,5 +1721,5 @@
 	IssmDouble fraction1,fraction2;
 	IssmDouble phi=element->GetGroundedPortion(xyz_list);
-	if(phi>0.00000001 && phi<0.99999999) partly_floating=true; 
+	if(phi>0.00000001 && phi<0.99999999) partly_floating=true;
 
 	int ig=-1;// needed for driving stress parameterization
@@ -1748,5 +1748,5 @@
 			/*Compute the discontinuous surface slope for this gauss point/subelement*/
 			this->ComputeSurfaceSlope(&slope[0],gauss_subelem,gauss,point1,fraction1,fraction2,ig,dim,element);
-		} 		
+		}
 		#endif
 
@@ -2130,5 +2130,5 @@
 			default:
 				_error_("index "<<point1<<" not supported yet");
-		}//}}}	
+		}//}}}
 	}
 	if(ig>3 && ig<8){ // GREEN element
@@ -2151,5 +2151,5 @@
 			default:
 				_error_("index "<<point1<<" not supported yet");
-		}//}}}	
+		}//}}}
 	}
 	if(ig>7){ // RED element
@@ -2172,5 +2172,5 @@
 			default:
 				_error_("index "<<point1<<" not supported yet");
-		}//}}}	
+		}//}}}
 	}
 
@@ -2198,5 +2198,5 @@
 }/*}}}*/
 void StressbalanceAnalysis::NodalFunctionsDerivativesRGB(IssmDouble* dbasis_subelem,Gauss* gauss_DG,Gauss* gauss_CG,int point1,IssmDouble fraction1,IssmDouble fraction2,int ig,int dim,Element* element){/*{{{*/
-	
+
 	/*Fetch number of nodes for this finite element*/
    int numnodes = element->GetNumberOfNodes();
@@ -2337,5 +2337,5 @@
    /*Get Jacobian*/
 	tria->GetJacobianInvert(&Jinv[0][0],xyz_list_RGB,gauss_CG); // gauss is not used here
-	
+
 	/*Build dbasis CG in the subelement*/
 	for(int i=0;i<numnodes;i++){
@@ -3435,5 +3435,5 @@
 	/*Get type of algorithm*/
 	int fe_FS;
-	bool isNitsche; 
+	bool isNitsche;
 
 	element->FindParam(&fe_FS,FlowequationFeFSEnum);
@@ -3760,5 +3760,5 @@
 	Input* vz_input = NULL;
 	if(dim==3){vz_input=element->GetInput(VzEnum); _assert_(vz_input);}
-	
+
 	/* prepare viscosity gradient for GLS */
 	IssmDouble NodalViscosity[6];
@@ -3914,5 +3914,5 @@
 	IssmDouble SU[3*(3+1)*6];
     Gauss* vert = element->NewGauss();
-	
+
 	element->FindParam(&FSreconditioning,StressbalanceFSreconditioningEnum);
 	/* Compute the nodal values of the viscosity */
@@ -3922,5 +3922,5 @@
 	}
 	delete vert;
-	
+
 	/* Start  looping on the number of gaussian points: */
 	Gauss* gauss=element->NewGauss(5);
@@ -3973,5 +3973,5 @@
 			pe->values[i*dim+0] += +rho_ice*forcex *Jdet*gauss->weight*vbasis[i];
 			pe->values[i*dim+1] += +rho_ice*forcey *Jdet*gauss->weight*vbasis[i];
-			if(dim==3) pe->values[i*dim+2] += +rho_ice*forcez*Jdet*gauss->weight*vbasis[i];	
+			if(dim==3) pe->values[i*dim+2] += +rho_ice*forcez*Jdet*gauss->weight*vbasis[i];
 
 			pe->values[i*dim+dim-1] += -rho_ice*gravity*Jdet*gauss->weight*vbasis[i];
@@ -4668,5 +4668,5 @@
 			pe->values[i*dim+0] += +rho_ice*forcex *Jdet*gauss->weight*vbasis[i];
 			pe->values[i*dim+1] += +rho_ice*forcey *Jdet*gauss->weight*vbasis[i];
-			if(dim==3) pe->values[i*dim+2] += +rho_ice*forcez*Jdet*gauss->weight*vbasis[i];	
+			if(dim==3) pe->values[i*dim+2] += +rho_ice*forcez*Jdet*gauss->weight*vbasis[i];
 
 			pe->values[i*dim+dim-1] += -rho_ice*gravity*Jdet*gauss->weight*vbasis[i];
@@ -4795,5 +4795,5 @@
 				}
 			}
-		}	
+		}
 		/* -------- Nitsche terms -------- */
 		/* n\cdot\nabla\phi*/
@@ -4814,5 +4814,5 @@
 			}
 		}
-		for(int i=0;i<pnumnodes;i++){	
+		for(int i=0;i<pnumnodes;i++){
 			nsigma[dim*vnumnodes+i] = -FSreconditioning*pbasis[i];
 		}
@@ -4832,5 +4832,5 @@
 				}
 			}
-		}	
+		}
 
 		/* pressure x velocity  component A12 */
Index: /issm/trunk-jpl/src/c/analyses/StressbalanceSIAAnalysis.cpp
===================================================================
--- /issm/trunk-jpl/src/c/analyses/StressbalanceSIAAnalysis.cpp	(revision 25553)
+++ /issm/trunk-jpl/src/c/analyses/StressbalanceSIAAnalysis.cpp	(revision 25554)
@@ -95,6 +95,6 @@
 
 	::CreateNodes(nodes,iomodel,StressbalanceSIAAnalysisEnum,P1Enum,SIAApproximationEnum);
-	for(int i=0;i<nodes->Size();i++){
-		Node* node=xDynamicCast<Node*>(nodes->GetObjectByOffset(i));
+	for(Object* & object: nodes->objects){
+		Node* node=xDynamicCast<Node*>(object);
 		int   sid = node->Sid();
 		int approximation=IoCodeToEnumVertexEquation(reCast<int>(iomodel->Data("md.flowequation.vertex_equation")[sid]));
@@ -602,5 +602,5 @@
 	for(i=0;i<numnodes;i++) vel[i]=sqrt(vx[i]*vx[i] + vy[i]*vy[i] + vz[i]*vz[i]);
 
-	/*For pressure: we have not computed pressure in this analysis, for this element. We are in 2D, 
+	/*For pressure: we have not computed pressure in this analysis, for this element. We are in 2D,
 	 *so the pressure is just the pressure at the bedrock: */
 	rho_ice  = element->FindParam(MaterialsRhoIceEnum);
@@ -612,5 +612,5 @@
 			for(i=0;i<numnodes;i++) pressure[i]=rho_ice*g*thickness[i];
 			break;
-		case Domain3DEnum:   
+		case Domain3DEnum:
 			element->GetVerticesCoordinates(&xyz_list);
 			element->GetInputListOnNodes(&surface[0],SurfaceEnum);
Index: /issm/trunk-jpl/src/c/classes/Constraints/Constraints.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/Constraints/Constraints.cpp	(revision 25553)
+++ /issm/trunk-jpl/src/c/classes/Constraints/Constraints.cpp	(revision 25554)
@@ -22,6 +22,6 @@
 void Constraints::ActivatePenaltyMethod(int in_analysis){/*{{{*/
 
-	for(int i=0;i<this->Size();i++){
-		Constraint* constraint=(Constraint*)this->GetObjectByOffset(i);
+	for(Object* & object: this->objects){
+		Constraint* constraint=(Constraint*)object;
 		constraint->ActivatePenaltyMethod();
 	}
Index: /issm/trunk-jpl/src/c/classes/Contours.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/Contours.cpp	(revision 25553)
+++ /issm/trunk-jpl/src/c/classes/Contours.cpp	(revision 25554)
@@ -37,8 +37,8 @@
 
 	/*open domain outline file for writing: */
-	if((fid=fopen(domainname,"w"))==NULL) _error_("could not open domain file " << domainname); 
-
-	for(int counter=0;counter<contours->Size();counter++){
-		contour=(Contour<double>*)contours->GetObjectByOffset(counter);
+	if((fid=fopen(domainname,"w"))==NULL) _error_("could not open domain file " << domainname);
+	int counter = 0;
+	for(Object* & object : contours->objects){
+		contour=(Contour<double>*)object;
 
 		/*Write header: */
@@ -56,4 +56,5 @@
 		/*Write blank line: */
 		if(counter<contours->Size()-1) fprintf(fid,"\n");
+		counter++;
 	}
 
Index: /issm/trunk-jpl/src/c/classes/Elements/Elements.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/Elements.cpp	(revision 25553)
+++ /issm/trunk-jpl/src/c/classes/Elements/Elements.cpp	(revision 25554)
@@ -51,7 +51,6 @@
 
 	/*Now go through all elements, and get how many nodes they own, unless they are clone nodes: */
-	for(int i=0;i<this->Size();i++){
-
-		Element* element=xDynamicCast<Element*>(this->GetObjectByOffset(i));
+	for(Object* & object : this->objects){
+		Element* element = xDynamicCast<Element*>(object);
 		numnodes=element->GetNumberOfNodes();
 		if(numnodes>max)max=numnodes;
Index: /issm/trunk-jpl/src/c/classes/ExternalResults/Results.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/ExternalResults/Results.cpp	(revision 25553)
+++ /issm/trunk-jpl/src/c/classes/ExternalResults/Results.cpp	(revision 25554)
@@ -33,6 +33,6 @@
 int Results::AddResult(ExternalResult* in_result){/*{{{*/
 
-	/*First, go through dataset of inputs and check whether any input 
-	 * with the same name is already in. If so, erase the corresponding 
+	/*First, go through dataset of inputs and check whether any input
+	 * with the same name is already in. If so, erase the corresponding
 	 * object before adding this new one: */
 
@@ -40,6 +40,6 @@
 	_assert_(in_result);
 
-	for(int i=0;i<this->Size();i++){
-		ExternalResult* result=xDynamicCast<ExternalResult*>(this->GetObjectByOffset(i));
+	for(Object* &object : this->objects){
+		ExternalResult* result=xDynamicCast<ExternalResult*>(object);
 
 		if(result->GetStep()==in_result->GetStep()){
@@ -64,7 +64,6 @@
 int Results::DeleteResult(int result_enum,int result_step){/*{{{*/
 
-	for(int i=0;i<this->Size();i++){
-		ExternalResult* result=xDynamicCast<ExternalResult*>(this->GetObjectByOffset(i));
-
+	for(Object* &object : this->objects){
+		ExternalResult* result=xDynamicCast<ExternalResult*>(object);
 		if(result->GetStep()==result_step){
 			if(strcmp(result->GetResultName(),EnumToStringx(result_enum))==0){
@@ -80,7 +79,6 @@
 ExternalResult* Results::FindResult(int result_enum){/*{{{*/
 
-	for(int i=0;i<this->Size();i++){
-		ExternalResult* result=xDynamicCast<ExternalResult*>(this->GetObjectByOffset(i));
-
+	for(Object* &object : this->objects){
+		ExternalResult* result=xDynamicCast<ExternalResult*>(object);
 		if(result->GetResultEnum()==result_enum){
 			return result;
@@ -99,6 +97,6 @@
 	parameters->FindParam(&io_gather,SettingsIoGatherEnum);
 
-	for(int i=0;i<this->Size();i++){
-		ExternalResult* result=xDynamicCast<ExternalResult*>(this->GetObjectByOffset(i));
+	for(Object* &object : this->objects){
+		ExternalResult* result=xDynamicCast<ExternalResult*>(object);
 		result->WriteData(fid,io_gather);
 	}
Index: /issm/trunk-jpl/src/c/classes/FemModel.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/FemModel.cpp	(revision 25553)
+++ /issm/trunk-jpl/src/c/classes/FemModel.cpp	(revision 25554)
@@ -108,5 +108,10 @@
 		this->parameters->FindParam(&isSSA,FlowequationIsSSAEnum);
 		this->parameters->FindParam(&domaintype,DomainTypeEnum);
-		for(int i=0;i<this->nummodels;i++) if(this->analysis_type_list[i]==StressbalanceAnalysisEnum){analysis_counter=i;break;}
+		for(int i=0;i<this->nummodels;i++) {
+			if(this->analysis_type_list[i]==StressbalanceAnalysisEnum){
+				analysis_counter=i;
+				break;
+			}
+		}
 		if(analysis_counter==-1) _error_("Could not find alias for analysis_type StressbalanceAnalysisEnum in list of FemModel analyses\n");
 		for(Object* & object : this->elements->objects){
@@ -1182,6 +1187,7 @@
 	IssmDouble* segmentlist    = xNewZeroInit<IssmDouble>(4*numseg);
 	IssmDouble* allsegmentlist = xNewZeroInit<IssmDouble>(4*numseg);
-	for(int i=0;i<segments->Size();i++){
-		Contour<IssmDouble>* segment=(Contour<IssmDouble>*)segments->GetObjectByOffset(i);
+	int i=0;
+	for(Object* & object : segments->objects){
+		Contour<IssmDouble>* segment=(Contour<IssmDouble>*)object;
 		_assert_(segment->nods == 2);
 		segmentlist[(numseg_offset+i)*4 + 0] = segment->x[0];
@@ -1189,4 +1195,5 @@
 		segmentlist[(numseg_offset+i)*4 + 2] = segment->x[1];
 		segmentlist[(numseg_offset+i)*4 + 3] = segment->y[1];
+		i++;
 	}
 
@@ -2295,6 +2302,7 @@
 
 			/*Go through our dependent variables, and compute the response:*/
-			for(int i=0;i<dependent_objects->Size();i++){
-				DependentObject* dep=(DependentObject*)dependent_objects->GetObjectByOffset(i);
+			int i = 0;
+			for(Object* & object : dependent_objects->objects){
+				DependentObject* dep=(DependentObject*)object;
 				dep->Responsex(&output_value,this);
 				if (my_rank==0) {
@@ -2307,4 +2315,5 @@
 					#endif
 				}
+				i++;
 			}
 		}
@@ -2437,7 +2446,7 @@
 
 							/*Fill-in vector*/
-							for(int j=0;j<this->loads->Size();j++){
-								if(this->loads->GetEnum(j)==ChannelEnum){
-									Channel* channel=(Channel*)this->loads->GetObjectByOffset(j);
+							for(Object* & object : this->loads->objects){
+								if(object->ObjectEnum()==ChannelEnum){
+									Channel* channel=(Channel*)object;
 									channel->WriteChannelCrossSection(values);
 								}
@@ -3200,6 +3209,6 @@
 	CreateNumberNodeToElementConnectivity(iomodel);
 	::CreateVertices(new_elements,new_vertices,iomodel,TransientSolutionEnum,true);
-	for(int i=0;i<new_vertices->Size();i++){
-		Vertex *vertex=(Vertex*)new_vertices->GetObjectByOffset(i);
+	for(Object* & object : new_vertices->objects){
+		Vertex *vertex=(Vertex*)object;
 		int     sid = vertex->Sid();
 		vertex->x=newx[sid];
@@ -3798,6 +3807,6 @@
 
 	/*Go through elements, and for each element, get vertices*/
-   for(int i=0;i<femmodel_elements->Size();i++){
-    	Element* element=xDynamicCast<Element*>(femmodel_elements->GetObjectByOffset(i));
+   for(Object* & object : femmodel_elements->objects){
+    	Element* element=xDynamicCast<Element*>(object);
     	element->GetVerticesSidList(elem_vertices);
     	vid1->SetValue(element->sid,elem_vertices[0],INS_VAL);
@@ -3972,8 +3981,6 @@
 
 	/*Get spcvx and spcvy of old mesh*/
-	for(int i=0;i<this->constraints_list[analysis_index]->Size();i++){
-
-		Constraint* constraint=(Constraint*)this->constraints_list[analysis_index]->GetObjectByOffset(i);
-
+	for(Object* & object : this->constraints_list[analysis_index]->objects){
+		Constraint* constraint=(Constraint*)object;
 		SpcStatic* spcstatic = xDynamicCast<SpcStatic*>(constraint);
 		int dof					= spcstatic->GetDof();
@@ -4890,11 +4897,11 @@
 	 * all DONE in Geographic coordinates: theta \in [-90,90], lambda \in [-180 180]
 	 */
-	for(int i=0;i<vertices->Size();i++){
-		int sid;
-		//Vertex* vertex=(Vertex*)vertices->GetObjectByOffset(i);
-		Vertex* vertex=xDynamicCast<Vertex*>(vertices->GetObjectByOffset(i));
-		sid=vertex->Sid();
-
-		lati=latitude[sid]/180*PI;	longi=longitude[sid]/180*PI; radi=radius[sid];
+	for(Object* & object : vertices->objects){
+		Vertex* vertex=xDynamicCast<Vertex*>(object);
+		int sid=vertex->Sid();
+
+		lati=latitude[sid]/180*PI;
+		longi=longitude[sid]/180*PI;
+		radi=radius[sid];
 
 		/*only first order terms are considered now: */
@@ -5315,9 +5322,7 @@
 		}
 		else{
-			int j=-1;
 			for(Object* & object : this->elements->objects){
 				/*Get the right transient input*/
 				Element* element = xDynamicCast<Element*>(object);
-				j+=1;
 				/*Get basal element*/
 				switch(domaintype){
@@ -5343,5 +5348,5 @@
 						element->GetInputListOnVertices(&values[0],input_enum[i]); //this is the enum to stack
 						element->GetVerticesLidList(vertexlids);
-						Penta* penta=xDynamicCast<Penta*>(elements->GetObjectByOffset(j));
+						Penta* penta=xDynamicCast<Penta*>(object);
 						for(;;){
 							transientinput->AddPentaTimeInput(subtime,numvertices,vertexlids,values,P1Enum);
Index: /issm/trunk-jpl/src/c/classes/Vertices.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/Vertices.cpp	(revision 25553)
+++ /issm/trunk-jpl/src/c/classes/Vertices.cpp	(revision 25554)
@@ -70,5 +70,5 @@
 	output->sorted_ids=NULL;
 	int objsize = this->numsorted;
-	if(this->sorted && objsize>0 && this->id_offsets){	
+	if(this->sorted && objsize>0 && this->id_offsets){
 		output->id_offsets=xNew<int>(objsize);
 		xMemCpy<int>(output->id_offsets,this->id_offsets,objsize);
@@ -154,6 +154,6 @@
 
 	/*Go through vertices, and for each vertex, object, report it cpu: */
-	for(int i=0;i<this->Size();i++){
-		Vertex* vertex=xDynamicCast<Vertex*>(this->GetObjectByOffset(i));
+	for(Object* & object : this->objects){
+      Vertex* vertex = xDynamicCast<Vertex*>(object);
 		lat->SetValue(vertex->sid,vertex->GetLatitude() ,INS_VAL);
 		lon->SetValue(vertex->sid,vertex->GetLongitude(),INS_VAL);
@@ -193,6 +193,6 @@
 	this->numberofvertices_local=this->Size();
 	this->numberofmasters_local=0;
-	for(int i=0;i<this->Size();i++){
-		Vertex* vertex=xDynamicCast<Vertex*>(this->GetObjectByOffset(i));
+	for(Object* & object : this->objects){
+      Vertex* vertex = xDynamicCast<Vertex*>(object);
 		if(!vertex->clone) this->numberofmasters_local++;
 	}
@@ -204,6 +204,6 @@
 
 	int lid = 0;
-	for(int i=0;i<this->Size();i++){
-		Vertex* vertex=xDynamicCast<Vertex*>(this->GetObjectByOffset(i));
+	for(Object* & object : this->objects){
+      Vertex* vertex = xDynamicCast<Vertex*>(object);
 		if(!vertex->clone){
 			vertex->lid=lid;
@@ -212,6 +212,6 @@
 		}
 	}
-	for(int i=0;i<this->Size();i++){
-		Vertex* vertex=xDynamicCast<Vertex*>(this->GetObjectByOffset(i));
+	for(Object* & object : this->objects){
+      Vertex* vertex = xDynamicCast<Vertex*>(object);
 		if(vertex->clone){
 			vertex->lid=lid;
@@ -229,6 +229,6 @@
 	xDelete<int>(all_num_masters);
 
-	for(int i=0;i<this->Size();i++){
-		Vertex* vertex=xDynamicCast<Vertex*>(this->GetObjectByOffset(i));
+	for(Object* & object : this->objects){
+      Vertex* vertex = xDynamicCast<Vertex*>(object);
 		vertex->pid = vertex->lid+offset;
 	}
Index: /issm/trunk-jpl/src/c/classes/kriging/Observations.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/kriging/Observations.cpp	(revision 25553)
+++ /issm/trunk-jpl/src/c/classes/kriging/Observations.cpp	(revision 25554)
@@ -283,5 +283,5 @@
 			}
 		}
-	}  
+	}
 
 	/*Assign output pointer*/
@@ -329,7 +329,9 @@
 		y   = xNew<IssmPDouble>(nobs);
 		obs = xNew<IssmPDouble>(nobs);
-		for(int i=0;i<this->Size();i++){
-			observation=xDynamicCast<Observation*>(this->GetObjectByOffset(i));
+		int i=0;
+		for(Object* & object: this->objects){
+			observation=xDynamicCast<Observation*>(object);
 			observation->WriteXYObs(&x[i],&y[i],&obs[i]);
+			i++;
 		}
 	}
@@ -455,5 +457,5 @@
 			if(stop) break;
 		}
-	}  
+	}
 	xDelete<IssmPDouble>(dists);
 	xDelete<int>(tempindices);
@@ -499,5 +501,5 @@
 	/*If we have less observations than mindata, return UNDEF*/
 	if(n_obs<mindata){
-		prediction = UNDEF; 
+		prediction = UNDEF;
 	}
 	else{
@@ -515,5 +517,5 @@
 			denominator += weight;
 		}
-		prediction = numerator/denominator; 
+		prediction = numerator/denominator;
 	}
 
@@ -543,6 +545,6 @@
 	/*If we have less observations than mindata, return UNDEF*/
 	if(n_obs<mindata){
-		*pprediction = -999.0; 
-		*perror      = -999.0; 
+		*pprediction = -999.0;
+		*perror      = -999.0;
 		return;
 	}
@@ -633,5 +635,5 @@
 	/*If we have less observations than mindata, return UNDEF*/
 	if(n_obs<mindata || n_obs<2){
-		prediction = UNDEF; 
+		prediction = UNDEF;
 	}
 	else{
@@ -714,9 +716,9 @@
 	for(j=0;j<n;j++) gamma[j]   = 0.0;
 
-	for(i=0;i<this->Size();i++){
-		observation1=xDynamicCast<Observation*>(this->GetObjectByOffset(i));
-
-		for(j=i+1;j<this->Size();j++){
-			observation2=xDynamicCast<Observation*>(this->GetObjectByOffset(j));
+	for(Object* & object : this->objects){
+		observation1=xDynamicCast<Observation*>(object);
+
+		for(Object* & object : this->objects){
+			observation2=xDynamicCast<Observation*>(object);
 
 			distance=sqrt(pow(observation1->x - observation2->x,2) + pow(observation1->y - observation2->y,2));
Index: /issm/trunk-jpl/src/c/cores/controladm1qn3_core.cpp
===================================================================
--- /issm/trunk-jpl/src/c/cores/controladm1qn3_core.cpp	(revision 25553)
+++ /issm/trunk-jpl/src/c/cores/controladm1qn3_core.cpp	(revision 25554)
@@ -1,5 +1,5 @@
 /*!\file: controladm1qn3_core.cpp
- * \brief: core of the control solution 
- */ 
+ * \brief: core of the control solution
+ */
 #include <ctime>
 #include <config.h>
@@ -281,6 +281,6 @@
 	codi_global.output_indices.clear();
 	#endif
-	for(int i=0;i<dependent_objects->Size();i++){
-		DependentObject* dep=(DependentObject*)dependent_objects->GetObjectByOffset(i);
+	for(Object* & object:dependent_objects->objects){
+		DependentObject* dep=(DependentObject*)object
 		if(solution_type==TransientSolutionEnum) output_value = dep->GetValue();
 		if(solution_type!=TransientSolutionEnum) dep->Responsex(&output_value,femmodel);
@@ -476,5 +476,5 @@
 	long         omode;
 	double		 f;
-	double		 dxmin,gttol; 
+	double		 dxmin,gttol;
 	int          maxsteps,maxiter;
 	int          intn,num_controls,num_cost_functions,solution_type;
@@ -596,5 +596,5 @@
 
 	for(int i=0;i<intn;i++) {
-		aX[i] = reCast<IssmDouble>(X[i]); 
+		aX[i] = reCast<IssmDouble>(X[i]);
 		aG[i] = reCast<IssmDouble>(G[i]);
 	}
Index: /issm/trunk-jpl/src/c/cores/controlvalidation_core.cpp
===================================================================
--- /issm/trunk-jpl/src/c/cores/controlvalidation_core.cpp	(revision 25553)
+++ /issm/trunk-jpl/src/c/cores/controlvalidation_core.cpp	(revision 25554)
@@ -1,5 +1,5 @@
 /*!\file: controlvalidation_core.cpp
- * \brief: core of the control solution 
- */ 
+ * \brief: core of the control solution
+ */
 
 #include "./cores.h"
@@ -217,6 +217,6 @@
 	dependents=xNew<IssmPDouble>(num_dependents);
 	codi_global.output_indices.clear();
-	for(int i=0;i<dependent_objects->Size();i++){
-		DependentObject* dep=(DependentObject*)dependent_objects->GetObjectByOffset(i);
+	for(Object* & object:dependent_objects->objects){
+		DependentObject* dep=(DependentObject*)object
 		if(solution_type==TransientSolutionEnum){
 			output_value = dep->GetValue();
@@ -235,5 +235,5 @@
 	j0 = J;
 	_printf0_("Initial cost function J(x) = "<<setw(12)<<setprecision(7)<<j0<<"\n");
-	_assert_(j0>0.); 
+	_assert_(j0>0.);
 	simul_stoptrace2();
 	/*initialize direction index in the weights vector: */
@@ -292,6 +292,6 @@
 		#if defined(_HAVE_CODIPACK_)
 		j=0.;
-		for(int i=0;i<dependent_objects->Size();i++){
-			DependentObject* dep=(DependentObject*)dependent_objects->GetObjectByOffset(i);
+		for(Object* & object:dependent_objects->objects){
+			DependentObject* dep=(DependentObject*)object
 			if(solution_type==TransientSolutionEnum){
 				output_value = dep->GetValue();
@@ -309,5 +309,5 @@
 		for(int i=0;i<n;i++) Den += alpha* G[i] * scaling_factors[0];
 		Ialpha = fabs((j - j0)/Den - 1.);
-		_assert_(fabs(Den)>0.); 
+		_assert_(fabs(Den)>0.);
 
 		_printf0_(" " << setw(11) << setprecision (5)<<alpha<<" " << setw(11) << setprecision (5)<<Ialpha<<"\n");
Index: /issm/trunk-jpl/src/c/cores/love_core.cpp
===================================================================
--- /issm/trunk-jpl/src/c/cores/love_core.cpp	(revision 25553)
+++ /issm/trunk-jpl/src/c/cores/love_core.cpp	(revision 25554)
@@ -1,5 +1,5 @@
 /*!\file: love_core.cpp
- * \brief: core of the LOVE numbers solution 
- */ 
+ * \brief: core of the LOVE numbers solution
+ */
 
 #include "./cores.h"
@@ -48,6 +48,6 @@
 	/*recover materials parameters: there is only one Matlitho, chase it down the hard way:*/
 	Matlitho* matlitho=NULL;
-	for (int i=femmodel->materials->Size()-1;i>=0;i--){
-		Material* material=xDynamicCast<Material*>(femmodel->materials->GetObjectByOffset(i));
+	for (Object* & object: femmodel->materials->objects){
+		Material* material=xDynamicCast<Material*>(object);
 		if (material->ObjectEnum()==MatlithoEnum)matlitho=xDynamicCast<Matlitho*>(material);
 	}
@@ -69,7 +69,7 @@
 	FourierLoveCorex(LoveKr,LoveKi,LoveHr,LoveHi,LoveLr,LoveLi,LoveKernelsReal,LoveKernelsImag,  //output
 			nfreq,frequencies,sh_nmax,sh_nmin,g0,r0,mu0,allow_layer_deletion,forcing_type,verbosemod, //parameter inputs
-			matlitho->numlayers, matlitho->radius, matlitho->viscosity, matlitho->lame_lambda, matlitho->lame_mu, 
+			matlitho->numlayers, matlitho->radius, matlitho->viscosity, matlitho->lame_lambda, matlitho->lame_mu,
 			matlitho->burgers_viscosity, matlitho->burgers_mu, matlitho->density, matlitho->isburgers, matlitho->issolid //matlitho inputs
-			); 
+			);
 
 	/*Add love matrices to results:*/
Index: /issm/trunk-jpl/src/c/cores/transient_core.cpp
===================================================================
--- /issm/trunk-jpl/src/c/cores/transient_core.cpp	(revision 25553)
+++ /issm/trunk-jpl/src/c/cores/transient_core.cpp	(revision 25554)
@@ -128,6 +128,6 @@
 		if(iscontrol && isautodiff){
 			/*Go through our dependent variables, and compute the response:*/
-			for(int i=0;i<dependent_objects->Size();i++){
-				DependentObject* dep=(DependentObject*)dependent_objects->GetObjectByOffset(i);
+			for(Object* & object:dependent_objects->objects){
+				DependentObject* dep=(DependentObject*)object;
 				IssmDouble  output_value;
 				dep->Responsex(&output_value,femmodel);
@@ -321,6 +321,6 @@
 
 		/*Go through our dependent variables, and compute the response:*/
-		for(int i=0;i<dependent_objects->Size();i++){
-			DependentObject* dep=(DependentObject*)dependent_objects->GetObjectByOffset(i);
+		for(Object* & object:dependent_objects->objects){
+			DependentObject* dep=(DependentObject*)object;
 			dep->Responsex(&output_value,femmodel);
 			dep->AddValue(output_value);
@@ -357,8 +357,9 @@
 
 	SetControlInputsFromVectorx(femmodel,X);
+
 	IssmDouble J = 0.;
 	if(0){
-		for(int i=0;i<dependent_objects->Size();i++){
-			DependentObject* dep=(DependentObject*)dependent_objects->GetObjectByOffset(i);
+		for(Object* & object:dependent_objects->objects){
+			DependentObject* dep=(DependentObject*)object;
 			J += dep->GetValue();
 		}
Index: /issm/trunk-jpl/src/c/main/kriging.cpp
===================================================================
--- /issm/trunk-jpl/src/c/main/kriging.cpp	(revision 25553)
+++ /issm/trunk-jpl/src/c/main/kriging.cpp	(revision 25554)
@@ -1,5 +1,5 @@
 /*!\file:  kriging.cpp
  * \brief: kriging main parallel program
- */ 
+ */
 
 #include "./issm.h"
@@ -57,6 +57,6 @@
 		results->AddObject(new GenericExternalResult<double*>(results->Size()+1,"predictions",predictions,ninterp,1,1,0));
 		results->AddObject(new GenericExternalResult<double*>(results->Size()+1,"error",error,ninterp,1,1,0));
-		for(int i=0;i<results->Size();i++){
-			ExternalResult* result=xDynamicCast<ExternalResult*>(results->GetObjectByOffset(i));
+		for(Object* & object :results->objects){
+			ExternalResult* result=xDynamicCast<ExternalResult*>(object);
 			result->WriteData(output_fid,1);
 		}
@@ -104,10 +104,10 @@
 
 	rootpatharg=argv[1];
-	if(strcmp(strstr(rootpatharg,"/"),"/")!=0){ 
+	if(strcmp(strstr(rootpatharg,"/"),"/")!=0){
 		rootpath       = xNew<char>(strlen(rootpatharg)+2); sprintf(rootpath,"%s/",rootpatharg);
-	}  
+	}
 	else{
 		rootpath       = xNew<char>(strlen(rootpatharg)+1); sprintf(rootpath,"%s",rootpatharg);
-	} 
+	}
 
 	modelname=argv[2];
Index: /issm/trunk-jpl/src/c/modules/AllocateSystemMatricesx/AllocateSystemMatricesx.cpp
===================================================================
--- /issm/trunk-jpl/src/c/modules/AllocateSystemMatricesx/AllocateSystemMatricesx.cpp	(revision 25553)
+++ /issm/trunk-jpl/src/c/modules/AllocateSystemMatricesx/AllocateSystemMatricesx.cpp	(revision 25554)
@@ -111,5 +111,5 @@
 	int localnumnodes       = femmodel->nodes->Size();
 	int numberofdofspernode = femmodel->nodes->MaxNumDofs(GsetEnum);
-	int M                   = femmodel->nodes->NumberOfDofs(set1enum);
+	//int M                   = femmodel->nodes->NumberOfDofs(set1enum);
 	int N                   = femmodel->nodes->NumberOfDofs(set2enum);
 	int m                   = femmodel->nodes->NumberOfDofsLocal(set1enum);
@@ -117,5 +117,4 @@
 	int numnodesperelement  = femmodel->elements->MaxNumNodes();
 	int numnodesperload     = femmodel->loads->MaxNumNodes();
-
 	int elementssize        = femmodel->elements->Size();
 	int loadssize           = femmodel->loads->Size();
Index: /issm/trunk-jpl/src/c/modules/ConfigureObjectsx/ConfigureObjectsx.cpp
===================================================================
--- /issm/trunk-jpl/src/c/modules/ConfigureObjectsx/ConfigureObjectsx.cpp	(revision 25553)
+++ /issm/trunk-jpl/src/c/modules/ConfigureObjectsx/ConfigureObjectsx.cpp	(revision 25554)
@@ -21,11 +21,11 @@
 	}
 	if(VerboseMProcessor()) _printf0_("      Configuring loads...\n");
-	for(int i=0;i<loads->Size();i++){
-		Load* load=(Load*)loads->GetObjectByOffset(i);
+	for(Object* object : loads->objects){
+		Load* load=(Load*)object;
 		load->Configure(elements,loads,nodes,vertices,materials,parameters);
 	}
 	if(VerboseMProcessor()) _printf0_("      Configuring materials...\n");
-	for(int i=0;i<materials->Size();i++){
-		Material* material=(Material*)materials->GetObjectByOffset(i);
+	for(Object* & object : materials->objects){
+		Material* material=(Material*)object;
 		material->Configure(elements);
 	}
Index: /issm/trunk-jpl/src/c/modules/ConstraintsStatex/ConstraintsStatex.cpp
===================================================================
--- /issm/trunk-jpl/src/c/modules/ConstraintsStatex/ConstraintsStatex.cpp	(revision 25553)
+++ /issm/trunk-jpl/src/c/modules/ConstraintsStatex/ConstraintsStatex.cpp	(revision 25554)
@@ -39,6 +39,6 @@
 	/*Deal with pengrid*/
 	if(femmodel->loads->numpenalties){
-		for(int i=0;i<femmodel->loads->Size();i++){
-			Load* load=(Load*)femmodel->loads->GetObjectByOffset(i);
+		for(Object* & object : femmodel->loads->objects){
+			Load* load=(Load*)object;
 			if(load->ObjectEnum()==PengridEnum){
 				Pengrid* pengrid=(Pengrid*)load;
@@ -50,5 +50,5 @@
 
 	ISSM_MPI_Reduce(&num_unstable_constraints,&sum_num_unstable_constraints,1,ISSM_MPI_INT,ISSM_MPI_SUM,0,IssmComm::GetComm() );
-	ISSM_MPI_Bcast(&sum_num_unstable_constraints,1,ISSM_MPI_INT,0,IssmComm::GetComm());                
+	ISSM_MPI_Bcast(&sum_num_unstable_constraints,1,ISSM_MPI_INT,0,IssmComm::GetComm());
 	num_unstable_constraints=sum_num_unstable_constraints;
 
Index: /issm/trunk-jpl/src/c/modules/ConstraintsStatex/RiftConstraintsState.cpp
===================================================================
--- /issm/trunk-jpl/src/c/modules/ConstraintsStatex/RiftConstraintsState.cpp	(revision 25553)
+++ /issm/trunk-jpl/src/c/modules/ConstraintsStatex/RiftConstraintsState.cpp	(revision 25554)
@@ -1,4 +1,4 @@
 /*!\file RiftConstraintsState.cpp
- * \brief: manage penalties for rifts 
+ * \brief: manage penalties for rifts
  */
 #include "./ConstraintsStateLocal.h"
@@ -32,6 +32,4 @@
 void RiftConstrain(int* pnum_unstable_constraints,Loads* loads,int configuration_type){
 
-	int			i;
-
 	/* generic object pointer: */
 	Riftfront* riftfront=NULL;
@@ -40,12 +38,10 @@
 	int unstable;
 	int sum_num_unstable_constraints;
-	int num_unstable_constraints=0;	
+	int num_unstable_constraints=0;
 
 	/*Enforce constraints: */
-	for (i=0;i<loads->Size();i++){
-
-		if (RiftfrontEnum==loads->GetEnum(i)){
-
-			load=(Load*)loads->GetObjectByOffset(i);
+	for (Object* & object : loads->objects){
+		if (RiftfrontEnum==object->ObjectEnum()){
+			load=(Load*)object;
 			riftfront=(Riftfront*)load;
 			riftfront->Constrain(&unstable);
@@ -55,5 +51,5 @@
 
 	ISSM_MPI_Reduce (&num_unstable_constraints,&sum_num_unstable_constraints,1,ISSM_MPI_INT,ISSM_MPI_SUM,0,IssmComm::GetComm() );
-	ISSM_MPI_Bcast(&sum_num_unstable_constraints,1,ISSM_MPI_INT,0,IssmComm::GetComm());                
+	ISSM_MPI_Bcast(&sum_num_unstable_constraints,1,ISSM_MPI_INT,0,IssmComm::GetComm());
 	num_unstable_constraints=sum_num_unstable_constraints;
 
@@ -75,9 +71,7 @@
 
 	/*Enforce constraints: */
-	for (i=0;i<loads->Size();i++){
-
-		if (RiftfrontEnum==loads->GetEnum(i)){
-
-			load=(Load*)loads->GetObjectByOffset(i);
+	for (Object* & object : loads->objects){
+		if (RiftfrontEnum==object->ObjectEnum()){
+			load=(Load*)object;
 			riftfront=(Riftfront*)load;
 			if (riftfront->IsFrozen()){
@@ -90,5 +84,5 @@
 	/*Is there just one found? that would mean we have frozen! : */
 	ISSM_MPI_Reduce (&found,&mpi_found,1,ISSM_MPI_INT,ISSM_MPI_MAX,0,IssmComm::GetComm() );
-	ISSM_MPI_Bcast(&mpi_found,1,ISSM_MPI_INT,0,IssmComm::GetComm());                
+	ISSM_MPI_Bcast(&mpi_found,1,ISSM_MPI_INT,0,IssmComm::GetComm());
 	found=mpi_found;
 
@@ -106,8 +100,7 @@
 
 	/*Enforce constraints: */
-	for (i=0;i<loads->Size();i++){
-
-		if (RiftfrontEnum==loads->GetEnum(i)){
-			load=(Load*)loads->GetObjectByOffset(i);
+	for (Object* & object : loads->objects){
+		if (RiftfrontEnum==object->ObjectEnum()){
+			load=(Load*)object;
 			riftfront=(Riftfront*)load;
 			riftfront->FreezeConstraints();
Index: /issm/trunk-jpl/src/c/modules/ContourToMeshx/ContourToMeshxt.cpp
===================================================================
--- /issm/trunk-jpl/src/c/modules/ContourToMeshx/ContourToMeshxt.cpp	(revision 25553)
+++ /issm/trunk-jpl/src/c/modules/ContourToMeshx/ContourToMeshxt.cpp	(revision 25554)
@@ -1,5 +1,5 @@
 /*!\file:  ContourToMeshxt.cpp
  * \brief  "thread" core code for interpolating values from a structured grid.
- */ 
+ */
 
 #ifdef HAVE_CONFIG_H
@@ -36,6 +36,6 @@
 
 	/*Loop through all contours: */
-	for (i=0;i<contours->Size();i++){
-		Contour<double>* contour=(Contour<double>*)contours->GetObjectByOffset(i);
+	for (Object* & object : contours->objects){
+		Contour<double>* contour=(Contour<double>*)object;
 		IsInPoly(in_nod,contour->x,contour->y,contour->nods,x,y,i0,i1,edgevalue);
 	}
Index: /issm/trunk-jpl/src/c/modules/ContourToNodesx/ContourToNodesx.cpp
===================================================================
--- /issm/trunk-jpl/src/c/modules/ContourToNodesx/ContourToNodesx.cpp	(revision 25553)
+++ /issm/trunk-jpl/src/c/modules/ContourToNodesx/ContourToNodesx.cpp	(revision 25554)
@@ -40,6 +40,6 @@
 	/*Loop through all contours: */
 	if(contours){
-		for(int i=0;i<contours->Size();i++){
-			Contour<IssmPDouble>* contour=(Contour<IssmPDouble>*)contours->GetObjectByOffset(i);
+		for(Object* & object:contours->objects){
+			Contour<IssmPDouble>* contour=(Contour<IssmPDouble>*)object;
 			IsInPoly(flags,contour->x,contour->y,contour->nods,x,y,0,nods,edgevalue);
 		}
Index: /issm/trunk-jpl/src/c/modules/CreateJacobianMatrixx/CreateJacobianMatrixx.cpp
===================================================================
--- /issm/trunk-jpl/src/c/modules/CreateJacobianMatrixx/CreateJacobianMatrixx.cpp	(revision 25553)
+++ /issm/trunk-jpl/src/c/modules/CreateJacobianMatrixx/CreateJacobianMatrixx.cpp	(revision 25554)
@@ -33,6 +33,6 @@
 		delete Je;
 	}
-	for (int i=0;i<femmodel->loads->Size();i++){
-		load=(Load*)femmodel->loads->GetObjectByOffset(i);
+	for (Object* & object : femmodel->loads->objects){
+		load=(Load*)object;
 		load->CreateJacobianMatrix(Jff);
 		load->PenaltyCreateJacobianMatrix(Jff,kmax);
Index: /issm/trunk-jpl/src/c/modules/CreateNodalConstraintsx/CreateNodalConstraintsx.cpp
===================================================================
--- /issm/trunk-jpl/src/c/modules/CreateNodalConstraintsx/CreateNodalConstraintsx.cpp	(revision 25553)
+++ /issm/trunk-jpl/src/c/modules/CreateNodalConstraintsx/CreateNodalConstraintsx.cpp	(revision 25554)
@@ -29,8 +29,8 @@
 	 ys=new Vector<IssmDouble>(slocalsize,ssize);
 
-	/*go through all nodes, and for the ones corresponding to this configuration_type, fill the 
+	/*go through all nodes, and for the ones corresponding to this configuration_type, fill the
 	 * constraints vector with the constraint values: */
-	for(int i=0;i<nodes->Size();i++){
-		Node* node=(Node*)nodes->GetObjectByOffset(i);
+	for(Object* & object: nodes->objects){
+		Node* node=xDynamicCast<Node*>(object);
 		node->CreateNodalConstraints(ys);
 	}
Index: /issm/trunk-jpl/src/c/modules/ExpToLevelSetx/ExpToLevelSetxt.cpp
===================================================================
--- /issm/trunk-jpl/src/c/modules/ExpToLevelSetx/ExpToLevelSetxt.cpp	(revision 25553)
+++ /issm/trunk-jpl/src/c/modules/ExpToLevelSetx/ExpToLevelSetxt.cpp	(revision 25554)
@@ -1,5 +1,5 @@
 /*!\file:  ExpToLevelSetxt.cpp
  * \brief  "thread" core code for figuring out level set value from a contour and a cloud of points.
- */ 
+ */
 
 #ifdef HAVE_CONFIG_H
@@ -38,6 +38,6 @@
 
 	/*Loop through all contours: */
-	for(i=0;i<contours->Size();i++){
-		Contour<double>* contour=(Contour<double>*)contours->GetObjectByOffset(i);
+	for(Object* & object : contours->objects){
+		Contour<double>* contour=(Contour<double>*)object;
 		ContourToLevelSet(distance,contour->x,contour->y,contour->nods,x,y,i0,i1);
 	}
@@ -73,9 +73,9 @@
 	 * We use the following notations:
 	 * segment: v=(x1,y1), w=(x2,y2)
-	 * point:   p=(x0,y0) 
+	 * point:   p=(x0,y0)
 	 */
 
    /*Get segment length square (avoid sqrt) |w-v|^2*/
-	double l2 = pow(x2-x1,2)+pow(y2-y1,2); 
+	double l2 = pow(x2-x1,2)+pow(y2-y1,2);
 
    /*segment is single point: v == w*/
Index: /issm/trunk-jpl/src/c/modules/KillIcebergsx/KillIcebergsx.cpp
===================================================================
--- /issm/trunk-jpl/src/c/modules/KillIcebergsx/KillIcebergsx.cpp	(revision 25553)
+++ /issm/trunk-jpl/src/c/modules/KillIcebergsx/KillIcebergsx.cpp	(revision 25554)
@@ -20,6 +20,7 @@
 
 	/*Step 1, go through all elements and put 1 in local_mask if the element is grounded*/
-	for(int i=0;i<femmodel->elements->Size();i++){
-		Element* element=xDynamicCast<Element*>(femmodel->elements->GetObjectByOffset(i));
+	int i=0;
+	for(Object* & object : femmodel->elements->objects){
+		Element* element=xDynamicCast<Element*>(object);
 		if(!element->IsIceInElement()){
 			/*Nothing to do, just flag element to speed up the computation*/
@@ -32,4 +33,5 @@
 			}
 		}
+		i++;
 	}
 
@@ -53,6 +55,7 @@
 			keepgoing    = false;
 			didsomething = 0;
-			for(int i=0;i<femmodel->elements->Size();i++){
-				Element* element=xDynamicCast<Element*>(femmodel->elements->GetObjectByOffset(i));
+			int i=0;
+			for(Object* & object : femmodel->elements->objects){
+				Element* element=xDynamicCast<Element*>(object);
 
 				if(!element_flag[i]){
@@ -78,4 +81,5 @@
 					}
 				}
+				i++;
 			}
 			iter++;
Index: /issm/trunk-jpl/src/c/modules/Mergesolutionfromftogx/Mergesolutionfromftogx.cpp
===================================================================
--- /issm/trunk-jpl/src/c/modules/Mergesolutionfromftogx/Mergesolutionfromftogx.cpp	(revision 25553)
+++ /issm/trunk-jpl/src/c/modules/Mergesolutionfromftogx/Mergesolutionfromftogx.cpp	(revision 25554)
@@ -17,5 +17,5 @@
 	int ssize=nodes->NumberOfDofs(SsetEnum);
 
-	/*serialize uf and ys: those two vectors will be indexed by the nodes, who are the only ones 
+	/*serialize uf and ys: those two vectors will be indexed by the nodes, who are the only ones
 	 *that know which values should be plugged into ug and where: */
 	if(ssize) if(flag_ys0) ys->Set(0.0);
@@ -33,6 +33,6 @@
 
 	/*Let nodes figure it out*/
-	for(int i=0;i<nodes->Size();i++){
-		Node* node=(Node*)nodes->GetObjectByOffset(i);
+	for(Object* & object: nodes->objects){
+		Node* node=xDynamicCast<Node*>(object);
 		node->VecMerge(ug,local_uf,indices_uf,local_ys,indices_ys);
 	}
Index: /issm/trunk-jpl/src/c/modules/OutputDefinitionsResponsex/OutputDefinitionsResponsex.cpp
===================================================================
--- /issm/trunk-jpl/src/c/modules/OutputDefinitionsResponsex/OutputDefinitionsResponsex.cpp	(revision 25553)
+++ /issm/trunk-jpl/src/c/modules/OutputDefinitionsResponsex/OutputDefinitionsResponsex.cpp	(revision 25554)
@@ -14,6 +14,6 @@
 
 	/*Now, go through the output definitions, and retrieve the object which corresponds to our requested response, output_string: */
-	for(int i=0;i<output_definitions->Size();i++){
-		Definition* definition=dynamic_cast<Definition*>(output_definitions->GetObjectByOffset(i));
+	for(Object* & object : output_definitions->objects){
+		Definition* definition=dynamic_cast<Definition*>(object);
 
 		char* name = definition->Name();
@@ -42,8 +42,6 @@
 
 	/*Now, go through the output definitions, and retrieve the object which corresponds to our requested response, output_enum: */
-	for(int i=0;i<output_definitions->Size();i++){
-
-		//Definition* definition=xDynamicCast<Definition*>(output_definitions->GetObjectByOffset(i));
-		Definition* definition=dynamic_cast<Definition*>(output_definitions->GetObjectByOffset(i));
+	for(Object* & object : output_definitions->objects){
+		Definition* definition=dynamic_cast<Definition*>(object);
 
 		int en = definition->DefinitionEnum();
@@ -59,5 +57,5 @@
 
 	/*If we are here, did not find the definition for this response, not good!: */
-	_error_("Could not find the response for output definition " << EnumToStringx(output_enum) 
+	_error_("Could not find the response for output definition " << EnumToStringx(output_enum)
 				<<" ("<<output_enum<<")"
 				<< " because could not find the definition itself!");
Index: /issm/trunk-jpl/src/c/modules/Reducevectorgtofx/Reducevectorgtofx.cpp
===================================================================
--- /issm/trunk-jpl/src/c/modules/Reducevectorgtofx/Reducevectorgtofx.cpp	(revision 25553)
+++ /issm/trunk-jpl/src/c/modules/Reducevectorgtofx/Reducevectorgtofx.cpp	(revision 25554)
@@ -1,4 +1,4 @@
 /*!\file Reducevectorgtofx
- * \brief reduce petsc vector from g set to s set (free dofs), using the nodeset partitioning 
+ * \brief reduce petsc vector from g set to s set (free dofs), using the nodeset partitioning
  * vectors.
  */
@@ -29,6 +29,6 @@
 
 	/*Let nodes figure it out*/
-	for(int i=0;i<nodes->Size();i++){
-		Node* node=(Node*)nodes->GetObjectByOffset(i);
+	for(Object* & object : nodes->objects){
+		Node* node=(Node*)object;
 		node->VecReduce(uf,local_ug,indices_ug);
 	}
Index: /issm/trunk-jpl/src/c/modules/ResetConstraintsx/ResetConstraintsx.cpp
===================================================================
--- /issm/trunk-jpl/src/c/modules/ResetConstraintsx/ResetConstraintsx.cpp	(revision 25553)
+++ /issm/trunk-jpl/src/c/modules/ResetConstraintsx/ResetConstraintsx.cpp	(revision 25554)
@@ -29,6 +29,6 @@
 
 	/*Reset pengrid to inactive mode*/
-	for(int i=0;i<femmodel->loads->Size();i++){
-		Load* load=(Load*)femmodel->loads->GetObjectByOffset(i);
+	for(Object* & object : femmodel->loads->objects){
+		Load* load=(Load*)object;
 		if(load->ObjectEnum()==PengridEnum){
 			Pengrid* pengrid=(Pengrid*)load;
@@ -48,6 +48,6 @@
 
 	/*Reset pengrid to inactive mode*/
-	for(int i=0;i<femmodel->loads->Size();i++){
-		Load* load=(Load*)femmodel->loads->GetObjectByOffset(i);
+	for(Object* & object: femmodel->loads->objects){
+		Load* load=(Load*)object;
 		if(load->ObjectEnum()==PengridEnum){
 			Pengrid* pengrid=(Pengrid*)load;
Index: /issm/trunk-jpl/src/c/modules/SpcNodesx/SpcNodesx.cpp
===================================================================
--- /issm/trunk-jpl/src/c/modules/SpcNodesx/SpcNodesx.cpp	(revision 25553)
+++ /issm/trunk-jpl/src/c/modules/SpcNodesx/SpcNodesx.cpp	(revision 25554)
@@ -10,6 +10,6 @@
 void SpcNodesx(Nodes* nodes,Constraints* constraints,Parameters* parameters){
 
-	for(int i=0;i<constraints->Size();i++){
-		Constraint* constraint=(Constraint*)constraints->GetObjectByOffset(i);
+	for(Object* & object: constraints->objects){
+		Constraint* constraint=xDynamicCast<Constraint*>(object);
 		constraint->ConstrainNode(nodes,parameters);
 	}
Index: /issm/trunk-jpl/src/c/modules/Trianglex/Trianglex.cpp
===================================================================
--- /issm/trunk-jpl/src/c/modules/Trianglex/Trianglex.cpp	(revision 25553)
+++ /issm/trunk-jpl/src/c/modules/Trianglex/Trianglex.cpp	(revision 25554)
@@ -42,10 +42,10 @@
 	/*Create initial triangulation to call triangulate(). First number of points:*/
 	in.numberofpoints=0;
-	for (i=0;i<domain->Size();i++){
-		contour=(Contour<IssmPDouble>*)domain->GetObjectByOffset(i);
+	for (Object* & object : domain->objects){
+		contour=(Contour<IssmPDouble>*)object;
 		in.numberofpoints+=contour->nods-1;
 	}
-	for (i=0;i<rifts->Size();i++){
-		contour=(Contour<IssmPDouble>*)rifts->GetObjectByOffset(i);
+	for (Object* & object : rifts->objects){
+		contour=(Contour<IssmPDouble>*)object;
 		in.numberofpoints+=contour->nods;
 	}
@@ -58,6 +58,6 @@
 
 	counter=0;
-	for (i=0;i<domain->Size();i++){
-		contour=(Contour<IssmPDouble>*)domain->GetObjectByOffset(i);
+	for (Object* & object : domain->objects){
+		contour=(Contour<IssmPDouble>*)object;
 		for (j=0;j<contour->nods-1;j++){
 			in.pointlist[2*counter+0]=contour->x[j];
@@ -66,6 +66,6 @@
 		}
 	}
-	for (i=0;i<rifts->Size();i++){
-		contour=(Contour<IssmPDouble>*)rifts->GetObjectByOffset(i);
+	for (Object* & object : rifts->objects){
+		contour=(Contour<IssmPDouble>*)object;
 		for (j=0;j<contour->nods;j++){
 			in.pointlist[2*counter+0]=contour->x[j];
@@ -85,10 +85,10 @@
 	/*Build segments. First figure out number of segments: holes and closed outlines have as many segments as vertices: */
 	in.numberofsegments=0;
-	for (i=0;i<domain->Size();i++){
-		contour=(Contour<IssmPDouble>*)domain->GetObjectByOffset(i);
+	for (Object* & object : domain->objects){
+		contour=(Contour<IssmPDouble>*)object;
 		in.numberofsegments+=contour->nods-1;
 	}
-	for(i=0;i<rifts->Size();i++){
-		contour=(Contour<IssmPDouble>*)rifts->GetObjectByOffset(i);
+	for (Object* & object : rifts->objects){
+		contour=(Contour<IssmPDouble>*)object;
 		/*for rifts, we have one less segment as we have vertices*/
 		in.numberofsegments+=contour->nods-1;
@@ -99,6 +99,6 @@
 	counter=0;
 	backcounter=0;
-	for (i=0;i<domain->Size();i++){
-		contour=(Contour<IssmPDouble>*)domain->GetObjectByOffset(i);
+	for (Object* & object : domain->objects){
+		contour=(Contour<IssmPDouble>*)object;
 		for (j=0;j<contour->nods-2;j++){
 			in.segmentlist[2*counter+0]=counter;
@@ -115,6 +115,6 @@
 	}
 	counter2=counter;
-	for (i=0;i<rifts->Size();i++){
-		contour=(Contour<IssmPDouble>*)rifts->GetObjectByOffset(i);
+	for (Object* & object : rifts->objects){
+		contour=(Contour<IssmPDouble>*)object;
 		for (j=0;j<(contour->nods-1);j++){
 			in.segmentlist[2*counter2+0]=counter;
Index: /issm/trunk-jpl/src/c/modules/UpdateDynamicConstraintsx/UpdateDynamicConstraintsx.cpp
===================================================================
--- /issm/trunk-jpl/src/c/modules/UpdateDynamicConstraintsx/UpdateDynamicConstraintsx.cpp	(revision 25553)
+++ /issm/trunk-jpl/src/c/modules/UpdateDynamicConstraintsx/UpdateDynamicConstraintsx.cpp	(revision 25554)
@@ -19,7 +19,7 @@
 	yg_serial=yg->ToMPISerial();
 
-	for(int i=0;i<constraints->Size();i++){
-		Constraint* constraint=(Constraint*)constraints->GetObjectByOffset(i);
-		if(constraint->ObjectEnum()==SpcDynamicEnum){ 
+	for(Object* & object : constraints->objects){
+		Constraint* constraint=(Constraint*)object;
+		if(constraint->ObjectEnum()==SpcDynamicEnum){
 			((SpcDynamic*)constraint)->SetDynamicConstraint(nodes,yg_serial);
 		}
Index: /issm/trunk-jpl/src/c/modules/VertexCoordinatesx/VertexCoordinatesx.cpp
===================================================================
--- /issm/trunk-jpl/src/c/modules/VertexCoordinatesx/VertexCoordinatesx.cpp	(revision 25553)
+++ /issm/trunk-jpl/src/c/modules/VertexCoordinatesx/VertexCoordinatesx.cpp	(revision 25554)
@@ -1,5 +1,5 @@
 /*!\file VertexCoordinatesx
- * \brief: compute a vector x,y and z of vertex coordinates by 
- * marching through all our vertices. 
+ * \brief: compute a vector x,y and z of vertex coordinates by
+ * marching through all our vertices.
  */
 
@@ -32,6 +32,6 @@
 
 	/*march through our vertices: */
-	for(i=0;i<vertices->Size();i++){
-		Vertex* vertex=(Vertex*)vertices->GetObjectByOffset(i);
+	for(Object* & object : vertices->objects){
+		Vertex* vertex=(Vertex*)object;
 		vertex->VertexCoordinates(vx,vy,vz,spherical);
 	}
Index: /issm/trunk-jpl/src/c/solutionsequences/solutionsequence_fct.cpp
===================================================================
--- /issm/trunk-jpl/src/c/solutionsequences/solutionsequence_fct.cpp	(revision 25553)
+++ /issm/trunk-jpl/src/c/solutionsequences/solutionsequence_fct.cpp	(revision 25554)
@@ -1,5 +1,5 @@
 /*!\file: solutionsequence_fct.cpp
  * \brief: numerical core of flux corrected transport solution
- */ 
+ */
 
 #include "../toolkits/toolkits.h"
@@ -109,6 +109,6 @@
 
 		dmax = dmax * 1.e+3;
-		for(int i=0;i<femmodel->constraints->Size();i++){
-			Constraint* constraint=(Constraint*)femmodel->constraints->GetObjectByOffset(i);
+		for(Object* & object: femmodel->constraints->objects){
+			Constraint* constraint=xDynamicCast<Constraint*>(object);
 			constraint->PenaltyDofAndValue(&dof,&d,femmodel->nodes,femmodel->parameters);
 			if(dof!=-1){
@@ -125,6 +125,6 @@
 
 		dmax = dmax * 1.e+3;
-		for(int i=0;i<femmodel->constraints->Size();i++){
-			Constraint* constraint=(Constraint*)femmodel->constraints->GetObjectByOffset(i);
+		for(Object* & object: femmodel->constraints->objects){
+			Constraint* constraint=xDynamicCast<Constraint*>(object);
 			constraint->PenaltyDofAndValue(&dof,&d,femmodel->nodes,femmodel->parameters);
 			if(dof!=-1){
@@ -180,6 +180,6 @@
 	if(USEPENALTYMETHOD){
 		/*Option 1: Penalty method*/
-		for(int i=0;i<femmodel->constraints->Size();i++){
-			Constraint* constraint=(Constraint*)femmodel->constraints->GetObjectByOffset(i);
+		for(Object* & object: femmodel->constraints->objects){
+			Constraint* constraint=xDynamicCast<Constraint*>(object);
 			constraint->PenaltyDofAndValue(&dof,&d,femmodel->nodes,femmodel->parameters);
 			d = d*dmax;
@@ -195,6 +195,6 @@
 		IssmDouble* rows_spc = xNew<IssmDouble>(numrows);
 		numrows = 0;
-		for(int i=0;i<femmodel->constraints->Size();i++){
-			Constraint* constraint=(Constraint*)femmodel->constraints->GetObjectByOffset(i);
+		for(Object* & object: femmodel->constraints->objects){
+			Constraint* constraint=xDynamicCast<Constraint*>(object);
 			constraint->PenaltyDofAndValue(&dof,&d,femmodel->nodes,femmodel->parameters);
 			if(dof!=-1){
@@ -459,5 +459,5 @@
 	/*Go solve lower order solution*/
 	femmodel->profiler->Start(SOLVER);
-	SolverxPetsc(&u,LHS,RHS,NULL,NULL, femmodel->parameters); 
+	SolverxPetsc(&u,LHS,RHS,NULL,NULL, femmodel->parameters);
 	femmodel->profiler->Stop(SOLVER);
 	MatFree(&LHS);
@@ -506,5 +506,5 @@
 
 	/*Update Element inputs*/
-	InputUpdateFromSolutionx(femmodel,uf); 
+	InputUpdateFromSolutionx(femmodel,uf);
 	delete uf;