Index: /issm/trunk/src/c/Container/Constraints.cpp
===================================================================
--- /issm/trunk/src/c/Container/Constraints.cpp	(revision 4235)
+++ /issm/trunk/src/c/Container/Constraints.cpp	(revision 4235)
@@ -0,0 +1,202 @@
+/*
+ * \file Constraints.c
+ * \brief: implementation of the Constraints class, derived from DataSet class
+ */
+
+/*Headers: {{{1*/
+#ifdef HAVE_CONFIG_H
+	#include "config.h"
+#else
+#error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!"
+#endif
+
+#include <vector>
+#include <functional>
+#include <algorithm>
+#include <iostream>
+
+#include "./DataSet.h"
+#include "../shared/shared.h"
+#include "../include/include.h"
+#include "../EnumDefinitions/EnumDefinitions.h"
+
+using namespace std;
+/*}}}*/
+
+/*Object constructors and destructor*/
+/*FUNCTION Constraints::Constraints(){{{1*/
+Constraints::Constraints(){
+	return;
+}
+/*}}}*/
+/*FUNCTION Constraints::Constraints(int in_enum){{{1*/
+Constraints::Constraints(int in_enum): DataSet(in_enum){
+	//do nothing;
+	return;
+}
+/*}}}*/
+/*FUNCTION Constraints::~Constraints(){{{1*/
+Constraints::~Constraints(){
+	return;
+}
+/*}}}*/
+
+/*Numerics: */
+/*FUNCTION Constraints::NumberOfLocalRgbs{{{1*/
+int   Constraints::NumberOfLocalRgbs(int analysis_type){
+
+	int i;
+	int  count=0;
+
+	for(i=0;i<this->Size();i++){
+
+		Object* object=(Object*)this->GetObjectByOffset(i);
+
+		/*Check this is a single point constraint (spc): */
+		if (object->Enum()==RgbEnum){
+
+			Rgb* rgb=(Rgb*)object;
+			if(rgb->InAnalysis(analysis_type))count++;
+		}
+	}
+
+	return count;
+}
+/*}}}*/
+/*FUNCTION Constraints::NumberOfConstraints{{{1*/
+int Constraints::NumberOfConstraints(void){
+
+	int localconstraints;
+	int numberofconstraints;
+
+	/*Get number of local constraints*/
+	localconstraints=this->Size();
+
+	/*figure out total number of constraints combining all the cpus (no clones here)*/
+	#ifdef _PARALLEL_
+	MPI_Reduce(&localconstraints,&numberofconstraints,1,MPI_INT,MPI_SUM,0,MPI_COMM_WORLD );
+	MPI_Bcast(&numberofconstraints,1,MPI_INT,0,MPI_COMM_WORLD);
+	#else
+	numberofconstraints=localconstraints;
+	#endif
+
+	return numberofconstraints;
+}
+/*}}}*/
+/*FUNCTION Constraints::SetupSpcs{{{1*/
+void   Constraints::SetupSpcs(Nodes* nodes,Vec yg,int analysis_type){
+
+	int i;
+
+	Node* node=NULL;
+	int nodeid;
+	int dof;
+	double value;
+
+	for(i=0;i<this->Size();i++){
+
+		Object* object=(Object*)this->GetObjectByOffset(i);
+
+		/*Check this is a single point constraint (spc): */
+		if(object->Enum()==SpcEnum){
+
+			Spc* spc=(Spc*)object;
+
+			if(spc->InAnalysis(analysis_type)){
+
+				/*Ok, this object is a constraint. Get the nodeid from the node it applies to: */
+				nodeid=spc->GetNodeId();
+				dof=spc->GetDof();
+				value=spc->GetValue();
+
+				/*Now, chase through nodes and find the corect node: */
+				node=(Node*)nodes->GetObjectById(NULL,nodeid);
+
+				/*Apply constraint: */
+				if(node){ //in case the spc is dealing with a node on another cpu
+					node->ApplyConstraint(yg,dof,value);
+				}
+			}
+
+		}
+	}
+
+	/*Assemble yg: */
+	VecAssemblyBegin(yg);
+	VecAssemblyEnd(yg);
+}
+/*}}}*/
+/*FUNCTION Constraints::SetupMpcs{{{1*/
+void Constraints::SetupMpcs(Mat Rmg,Nodes* nodes,int analysis_type){
+
+	int i;
+
+	int  nodeid1;
+	int  nodeid2;
+	int  dof;
+
+	int  dof1;
+	int  dof2;
+
+
+	Node* node1=NULL;
+	Node* node2=NULL;
+
+	int count=-1;
+
+	for(i=0;i<this->Size();i++){
+
+		Object* object=(Object*)this->GetObjectByOffset(i);
+
+		/*Check this is a mutiple point constraint (spc): */
+		if(object->Enum()==RgbEnum){
+
+			Rgb* rgb=(Rgb*)object;
+
+			if(rgb->InAnalysis(analysis_type)){
+
+				/*we found an rgb, increment counter, so that row index for Rmg is up to date: */
+				count++;
+
+
+				nodeid1=rgb->GetNodeId1();
+				nodeid2=rgb->GetNodeId2();
+				dof=rgb->GetDof();
+
+				/*For this rgb, find the nodes that go with it: */
+				node1=(Node*)nodes->GetObjectById(NULL,nodeid1);
+				node2=(Node*)nodes->GetObjectById(NULL,nodeid2);
+
+				if ((node1 && !node2) || (!node1 && node2)){
+					/*we are missing one node, not good!*/
+					ISSMERROR("%s%p%s%p"," in Rgb, missing one node. node1: ",node1," node2: ",node2);
+				}
+
+				if(!node1 && !node2){
+					/*That's ok, this Rgb can't find those nodes, so leave them alone. They are probably not on this 
+					 * cpu!*/
+				}
+				else{
+					/*Ok, this cpu owns both nodes. Put dof for node1 into m set, unless it is already there, 
+					 * in which case node2 gets into the m set: */
+					if(node1->DofIsInMSet(dof-1)){
+						node2->DofInMSet(dof-1);
+					}
+					else{
+						node1->DofInMSet(dof-1);
+					}
+
+					/*Plug values into Rmg. We essentially want dofs from node1 and node2 to be the 
+					 *same: */
+					dof1=node1->GetDof(dof-1); //matlab indexing
+					dof2=node2->GetDof(dof-1); //matlab indexing
+
+					MatSetValue(Rmg,count,dof1,1.0,INSERT_VALUES);
+					MatSetValue(Rmg,count,dof2,-1.0,INSERT_VALUES);
+
+				}
+			}
+		}
+	}
+}
+/*}}}*/
Index: /issm/trunk/src/c/Container/DataSet.cpp
===================================================================
--- /issm/trunk/src/c/Container/DataSet.cpp	(revision 4235)
+++ /issm/trunk/src/c/Container/DataSet.cpp	(revision 4235)
@@ -0,0 +1,561 @@
+/*
+ * \file DataSet.c
+ * \brief: implementation of the DataSet class, and derived classes Inputs and Parameters
+ */
+
+/*Headers: {{{1*/
+#ifdef HAVE_CONFIG_H
+	#include "config.h"
+#else
+#error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!"
+#endif
+
+#include <vector>
+#include <functional>
+#include <algorithm>
+#include <iostream>
+
+#include "./DataSet.h"
+#include "../objects/objects.h"
+#include "../shared/shared.h"
+#include "../include/include.h"
+#include "../EnumDefinitions/EnumDefinitions.h"
+
+using namespace std;
+/*}}}*/
+
+/*Constructors/Destructors*/
+/*FUNCTION DataSet::DataSet(){{{1*/
+DataSet::DataSet(){
+	
+	sorted=0;
+	sorted_ids=NULL;
+	id_offsets=NULL;
+
+}
+/*}}}*/
+/*FUNCTION DataSet::DataSet(int dataset_enum){{{1*/
+DataSet::DataSet(int dataset_enum){
+	enum_type=dataset_enum;
+	
+	sorted=0;
+	sorted_ids=NULL;
+	id_offsets=NULL;
+
+}
+/*}}}*/
+/*FUNCTION DataSet::Copy{{{1*/
+DataSet*   DataSet::Copy(void){
+
+
+	DataSet* copy=NULL;
+	vector<Object*>::iterator object;
+	Object* object_copy=NULL;
+
+	copy=new DataSet(enum_type);
+
+	copy->sorted=sorted;
+	copy->presorted=presorted;
+	if(sorted_ids){
+		copy->sorted_ids=(int*)xmalloc(objects.size()*sizeof(int));
+		memcpy(copy->sorted_ids,sorted_ids,objects.size()*sizeof(int));
+	}
+	if(id_offsets){
+		copy->id_offsets=(int*)xmalloc(objects.size()*sizeof(int));
+		memcpy(copy->id_offsets,id_offsets,objects.size()*sizeof(int));
+	}
+
+	/*Now we need to deep copy the objects: */
+	for ( object=objects.begin() ; object < objects.end(); object++ ){
+
+		/*Call echo on object: */
+		object_copy = (*object)->copy();
+		copy->AddObject(object_copy);
+	}
+	return copy;
+}
+/*}}}*/
+/*FUNCTION DataSet::~DataSet{{{1*/
+DataSet::~DataSet(){
+	clear();
+	xfree((void**)&sorted_ids);
+	xfree((void**)&id_offsets);
+}
+/*}}}*/
+
+/*I/O*/
+/*FUNCTION DataSet::Marshall{{{1*/
+char* DataSet::Marshall(){
+
+	vector<Object*>::iterator object;
+	int                       object_size;
+	int                       marshalled_dataset_size=0;
+	char*                     marshalled_dataset=NULL;
+	char*                     old_marshalled_dataset=NULL;
+
+	/*First get size of marshalled dataset: */
+	object_size=(int)objects.size();
+
+	marshalled_dataset_size=MarshallSize();
+	
+	/*Allocate marshalled dataset: */
+	marshalled_dataset=(char*)xmalloc(marshalled_dataset_size*sizeof(char)); 
+
+	/*Keep track of old_marshalled_dataset: */
+	old_marshalled_dataset=marshalled_dataset;
+
+	/*Store internals of dataset first: */
+	memcpy(marshalled_dataset,&object_size,sizeof(int)); marshalled_dataset+=sizeof(int);
+	memcpy(marshalled_dataset,&sorted,sizeof(int)); marshalled_dataset+=sizeof(int);
+	if(sorted){
+		if(object_size)memcpy(marshalled_dataset,sorted_ids,object_size*sizeof(int)); marshalled_dataset+=object_size*sizeof(int);
+		if(object_size)memcpy(marshalled_dataset,id_offsets,object_size*sizeof(int)); marshalled_dataset+=object_size*sizeof(int);
+	}
+
+	for ( object=objects.begin() ; object < objects.end(); object++ ){
+		(*object)->Marshall(&marshalled_dataset);
+	}
+
+	/* Ok, marshalled_dataset now points to the end of the original marshalled_dataset pointer 
+	 * before  we started the loop on objects. Get object to point right again: */
+	marshalled_dataset-=marshalled_dataset_size;
+
+	/*We should be back to old_marshalled_dataset: check and abort if that's not the case, 
+	 * because this is a nasty error: */
+	if (marshalled_dataset!=old_marshalled_dataset){
+		ISSMERROR("final marshalled dataset is different from initial one!"); 
+		abort();
+	}
+
+	/*Return: */
+	return marshalled_dataset;
+}
+/*}}}*/
+/*FUNCTION DataSet::MarshallSize{{{1*/
+int DataSet::MarshallSize(){
+
+	vector<Object*>::iterator object;
+	int                      marshalled_dataset_size=0;
+
+
+	for ( object=objects.begin() ; object < objects.end(); object++ ){
+		marshalled_dataset_size+= (*object)->MarshallSize();
+	}
+
+	marshalled_dataset_size+=sizeof(int); //objects size
+	marshalled_dataset_size+=sizeof(int); //sorted size
+	if(sorted){
+		marshalled_dataset_size+=(int)objects.size()*sizeof(int); //sorted ids
+		marshalled_dataset_size+=(int)objects.size()*sizeof(int); //id offsets
+	}
+
+	return marshalled_dataset_size;
+}
+/*}}}*/
+/*FUNCTION DataSet::Demarshall{{{1*/
+DataSet* DataSetDemarshall(char* marshalled_dataset){
+
+	return DataSetDemarshallRaw(&marshalled_dataset);
+
+}
+/*}}}*/
+/*FUNCTION DataSet::DemarshallRaw{{{1*/
+DataSet* DataSetDemarshallRaw(char** pmarshalled_dataset){
+
+	int i;
+
+	DataSet* dataset=NULL;
+	int      numobjects=0;
+	int      enum_type;
+	Object*  object=NULL;
+	int      sorted;
+	int*     sorted_ids=NULL;
+	int*     id_offsets=NULL;
+	char*    marshalled_dataset=NULL;
+
+	/*recover marshalled_dataset pointer: */
+	marshalled_dataset=*pmarshalled_dataset;
+
+	/*initialize dataset: */
+	dataset=new DataSet();
+
+	/*Get internals first: */
+	memcpy(&numobjects,marshalled_dataset,sizeof(int)); marshalled_dataset+=sizeof(int);
+	memcpy(&sorted,marshalled_dataset,sizeof(int)); marshalled_dataset+=sizeof(int);
+	if(sorted){
+		if(numobjects){
+			sorted_ids=(int*)xmalloc(numobjects*sizeof(int));
+			id_offsets=(int*)xmalloc(numobjects*sizeof(int));
+			memcpy(sorted_ids,marshalled_dataset,numobjects*sizeof(int)); marshalled_dataset+=numobjects*sizeof(int);
+			memcpy(id_offsets,marshalled_dataset,numobjects*sizeof(int)); marshalled_dataset+=numobjects*sizeof(int);
+		}
+		dataset->SetSorting(sorted_ids,id_offsets);
+	}
+
+	for(i=0;i<numobjects;i++){
+
+		/*get enum type of object: */
+		memcpy(&enum_type,marshalled_dataset,sizeof(int)); marshalled_dataset+=sizeof(int);
+
+		if(enum_type==NodeEnum){
+			Node* node=NULL;
+			node=new Node();
+			node->Demarshall(&marshalled_dataset);
+			dataset->AddObject(node);
+		}
+		else if(enum_type==VertexEnum){
+			Vertex* vertex=NULL;
+			vertex=new Vertex();
+			vertex->Demarshall(&marshalled_dataset);
+			dataset->AddObject(vertex);
+		}
+		else if(enum_type==DoubleParamEnum){
+			DoubleParam* doubleparam=NULL;
+			doubleparam=new DoubleParam();
+			doubleparam->Demarshall(&marshalled_dataset);
+			dataset->AddObject(doubleparam);
+		}
+		else if(enum_type==TriaEnum){
+			Tria* tria=NULL;
+			tria=new Tria();
+			tria->Demarshall(&marshalled_dataset);
+			dataset->AddObject(tria);
+		}
+		else if(enum_type==TriaVertexInputEnum){
+			TriaVertexInput* triavertexinput=NULL;
+			triavertexinput=new TriaVertexInput();
+			triavertexinput->Demarshall(&marshalled_dataset);
+			dataset->AddObject(triavertexinput);
+		}
+		else if(enum_type==PentaVertexInputEnum){
+			PentaVertexInput* pentavertexinput=NULL;
+			pentavertexinput=new PentaVertexInput();
+			pentavertexinput->Demarshall(&marshalled_dataset);
+			dataset->AddObject(pentavertexinput);
+		}
+		else if(enum_type==SingVertexInputEnum){
+			SingVertexInput* singvertexinput=NULL;
+			singvertexinput=new SingVertexInput();
+			singvertexinput->Demarshall(&marshalled_dataset);
+			dataset->AddObject(singvertexinput);
+		}
+		else if(enum_type==BeamVertexInputEnum){
+			BeamVertexInput* beamvertexinput=NULL;
+			beamvertexinput=new BeamVertexInput();
+			beamvertexinput->Demarshall(&marshalled_dataset);
+			dataset->AddObject(beamvertexinput);
+		}
+		else if(enum_type==SingEnum){
+			Sing* sing=NULL;
+			sing=new Sing();
+			sing->Demarshall(&marshalled_dataset);
+			dataset->AddObject(sing);
+		}
+		else if(enum_type==BeamEnum){
+			Beam* beam=NULL;
+			beam=new Beam();
+			beam->Demarshall(&marshalled_dataset);
+			dataset->AddObject(beam);
+		}
+		else if(enum_type==PentaEnum){
+			Penta* penta=NULL;
+			penta=new Penta();
+			penta->Demarshall(&marshalled_dataset);
+			dataset->AddObject(penta);
+		}
+		else if(enum_type==MaticeEnum){
+			Matice* matice=NULL;
+			matice=new Matice();
+			matice->Demarshall(&marshalled_dataset);
+			dataset->AddObject(matice);
+		}
+		else if(enum_type==MatparEnum){
+			Matpar* matpar=NULL;
+			matpar=new Matpar();
+			matpar->Demarshall(&marshalled_dataset);
+			dataset->AddObject(matpar);
+		}
+		else if(enum_type==SpcEnum){
+			Spc* spc=NULL;
+			spc=new Spc();
+			spc->Demarshall(&marshalled_dataset);
+			dataset->AddObject(spc);
+		}
+		else if(enum_type==PengridEnum){
+			Pengrid* pengrid=NULL;
+			pengrid=new Pengrid();
+			pengrid->Demarshall(&marshalled_dataset);
+			dataset->AddObject(pengrid);
+		}
+		else if(enum_type==PenpairEnum){
+			Penpair* penpair=NULL;
+			penpair=new Penpair();
+			penpair->Demarshall(&marshalled_dataset);
+			dataset->AddObject(penpair);
+		}
+		else if(enum_type==IcefrontEnum){
+			Icefront* icefront=NULL;
+			icefront=new Icefront();
+			icefront->Demarshall(&marshalled_dataset);
+			dataset->AddObject(icefront);
+		}
+		else if(enum_type==NumericalfluxEnum){
+			Numericalflux* numericalflux=NULL;
+			numericalflux=new Numericalflux();
+			numericalflux->Demarshall(&marshalled_dataset);
+			dataset->AddObject(numericalflux);
+		}
+		else if(enum_type==RgbEnum){
+			Rgb* rgb=NULL;
+			rgb=new Rgb();
+			rgb->Demarshall(&marshalled_dataset);
+			dataset->AddObject(rgb);
+		}
+		else if(enum_type==RiftfrontEnum){
+			Riftfront* riftfront=NULL;
+			riftfront=new Riftfront();
+			riftfront->Demarshall(&marshalled_dataset);
+			dataset->AddObject(riftfront);
+		}
+		else if(enum_type==DoubleInputEnum){
+			DoubleInput* doubleinput=NULL;
+			doubleinput=new DoubleInput();
+			doubleinput->Demarshall(&marshalled_dataset);
+			dataset->AddObject(doubleinput);
+		}
+		else if(enum_type==IntInputEnum){
+			IntInput* intinput=NULL;
+			intinput=new IntInput();
+			intinput->Demarshall(&marshalled_dataset);
+			dataset->AddObject(intinput);
+		}
+		else if(enum_type==BoolInputEnum){
+			BoolInput* boolinput=NULL;
+			boolinput=new BoolInput();
+			boolinput->Demarshall(&marshalled_dataset);
+			dataset->AddObject(boolinput);
+		}
+		else if(enum_type==IntParamEnum){
+			IntParam* intparam=NULL;
+			intparam=new IntParam();
+			intparam->Demarshall(&marshalled_dataset);
+			dataset->AddObject(intparam);
+		}
+		else if(enum_type==BoolParamEnum){
+			BoolParam* boolparam=NULL;
+			boolparam=new BoolParam();
+			boolparam->Demarshall(&marshalled_dataset);
+			dataset->AddObject(boolparam);
+		}
+		else if(enum_type==StringParamEnum){
+			StringParam* stringparam=NULL;
+			stringparam=new StringParam();
+			stringparam->Demarshall(&marshalled_dataset);
+			dataset->AddObject(stringparam);
+		}
+		else{
+			ISSMERROR("could not recognize enum type: %i (%s)",enum_type,EnumAsString(enum_type));
+		}
+
+	}
+
+	/*Assign output pointers:*/
+	*pmarshalled_dataset=marshalled_dataset;
+	
+	return dataset;
+}
+/*}}}*/
+
+/*Specific methods*/
+/*FUNCTION DataSet::AddObject{{{1*/
+int  DataSet::AddObject(Object* object){
+
+	objects.push_back(object);
+
+	return 1;
+}
+/*}}}*/
+/*FUNCTION DataSet::clear{{{1*/
+void  DataSet::clear(){
+
+	vector<Object*>::iterator object;
+
+	for ( object=objects.begin() ; object < objects.end(); object++ ){
+		delete (*object);
+	}
+	objects.clear();
+}
+/*}}}*/
+/*FUNCTION DataSet::DeleteObject{{{1*/
+int  DataSet::DeleteObject(Object* object){
+
+	vector<Object*>::iterator iterator;
+
+	if(object){
+		iterator = find(objects.begin(), objects.end(),object);
+		delete *iterator;
+		objects.erase(iterator);
+	}
+
+}
+/*}}}*/
+/*FUNCTION DataSet::DeepEcho{{{1*/
+void DataSet::DeepEcho(){
+
+
+	vector<Object*>::iterator object;
+
+	if(this==NULL)ISSMERROR(" trying to echo a NULL dataset");
+
+	_printf_("DataSet echo: %i objects\n",objects.size());
+
+	for ( object=objects.begin() ; object < objects.end(); object++ ){
+
+		/*Call deep echo on object: */
+		(*object)->DeepEcho();
+
+	}
+	return;
+}
+/*}}}*/
+/*FUNCTION DataSet::Echo{{{1*/
+void DataSet::Echo(){
+
+
+	vector<Object*>::iterator object;
+
+	if(this==NULL)ISSMERROR(" trying to echo a NULL dataset");
+
+	_printf_("DataSet echo: %i objects\n",objects.size());
+
+	for ( object=objects.begin() ; object < objects.end(); object++ ){
+
+		/*Call echo on object: */
+		(*object)->Echo();
+
+	}
+	return;
+}
+/*}}}*/
+/*FUNCTION DataSet::GetEnum{{{1*/
+int  DataSet::GetEnum(){
+	return enum_type;
+}
+/*}}}*/
+/*FUNCTION DataSet::GetEnum(int offset){{{1*/
+int   DataSet::GetEnum(int offset){
+
+	return objects[offset]->Enum();
+
+}
+/*}}}*/
+/*FUNCTION DataSet::GetObjectByOffset{{{1*/
+Object* DataSet::GetObjectByOffset(int offset){
+
+	return objects[offset];
+
+}
+/*}}}*/
+/*FUNCTION DataSet::GetObjectById{{{1*/
+Object* DataSet::GetObjectById(int* poffset,int eid){
+
+	int id_offset;
+	int offset;
+	int i;
+
+	if(!sorted)ISSMERROR(" trying to binary search on a non-sorted dataset!");
+
+	/*Carry out a binary search on the sorted_ids: */
+	if(!binary_search(&id_offset,eid, sorted_ids,objects.size())){
+		ISSMERROR("could not find object with id %i in DataSet %s",eid,EnumAsString(enum_type));
+	}
+
+	/*Convert  the id offset into sorted offset: */
+	offset=id_offsets[id_offset];
+
+	/*Assign output pointers if requested:*/
+	if (poffset)*poffset=offset;
+
+	/*Return object at offset position in objects :*/
+	return objects[offset];
+}
+/*}}}*/
+/*FUNCTION DataSet::Presort{{{1*/
+void DataSet::Presort(){
+
+	/*vector of objects is already sorted, just allocate the sorted ids and their
+	 * offsets:*/
+	int i;
+
+	if(objects.size()){
+
+		/*Delete existing ids*/
+		xfree((void**)&sorted_ids);
+		xfree((void**)&id_offsets);
+
+		/*Allocate new ids*/
+		sorted_ids=(int*)xmalloc(objects.size()*sizeof(int));
+		id_offsets=(int*)xmalloc(objects.size()*sizeof(int));
+
+		/*Build id_offsets and sorted_ids*/
+		for(i=0;i<objects.size();i++){
+			id_offsets[i]=i;
+			sorted_ids[i]=objects[i]->Id();
+		}
+	}
+
+	/*set sorted flag: */
+	sorted=1;
+}
+/*}}}*/
+/*FUNCTION DataSet::SetSorting{{{1*/
+void DataSet::SetSorting(int* in_sorted_ids,int* in_id_offsets){
+
+	sorted=1;
+	sorted_ids=in_sorted_ids;
+	id_offsets=in_id_offsets;
+}
+/*}}}*/
+/*FUNCTION DataSet::Size{{{1*/
+int  DataSet::Size(void){
+
+	return objects.size();
+}
+/*}}}*/
+/*FUNCTION DataSet::Sort{{{1*/
+void DataSet::Sort(){
+
+	/*Only sort if we are not already sorted: */
+	if(!sorted){
+		ISSMERROR(" not implemented yet!");
+	}
+}
+/*}}}*/
+/*FUNCTION DataSet::Configure{{{1*/
+void DataSet::Configure(Elements* elements,Loads* loads, DataSet* nodes, Vertices* vertices, Materials* materials,Parameters* parameters){
+
+	vector<Object*>::iterator object;
+	Element* element=NULL;
+	Load* load=NULL;
+	Node* node=NULL;
+
+	for ( object=objects.begin() ; object < objects.end(); object++ ){
+
+		if(EnumIsElement((*object)->Enum())){
+
+			element=(Element*)(*object);
+			element->Configure(elements,loads,nodes,materials,parameters);
+		}
+		if(EnumIsLoad((*object)->Enum())){
+			load=(Load*)(*object);
+			load->Configure(elements,loads,nodes,vertices,materials,parameters);
+		}
+
+		if((*object)->Enum()==NodeEnum){
+			node=(Node*)(*object);
+			node->Configure(nodes,vertices);
+		}
+	}
+
+}
+/*}}}*/
Index: /issm/trunk/src/c/Container/DataSet.h
===================================================================
--- /issm/trunk/src/c/Container/DataSet.h	(revision 4235)
+++ /issm/trunk/src/c/Container/DataSet.h	(revision 4235)
@@ -0,0 +1,298 @@
+/*
+ * DataSet.h: declaration of DataSet,Parameters and Inputs classes
+ */
+
+#ifndef _DATASET_H_
+#define _DATASET_H_
+
+#include <vector>
+#include "../objects/Object.h"
+#include "../toolkits/toolkits.h"
+#include "../EnumDefinitions/EnumDefinitions.h"
+
+/*forward declarations */
+class Constraints;
+class Materials;
+class Parameters;
+class Elements;
+class Vertices;
+class Loads;
+class Nodes;
+class DataSet;
+class Inputs;
+
+
+/********************************************************DATASET************************************************/
+class DataSet{
+	
+	public: 
+		
+		/*internals: */
+		std::vector<Object*> objects;
+		
+		/*type of dataset: */
+		int             enum_type;
+		
+		/*sorting: */
+		int             sorted;
+		int             presorted;
+		int*            sorted_ids;
+		int*            id_offsets;
+
+		/*constructors, destructors: {{{1*/
+		DataSet();
+		DataSet(int enum_type);
+		~DataSet();
+		/*}}}*/
+		/*management: {{{1*/
+		int   GetEnum();
+		int   GetEnum(int offset);
+		void  Echo();
+		void  DeepEcho();
+		char* Marshall();
+		int   MarshallSize();
+		int   AddObject(Object* object);
+		int   DeleteObject(int id);
+		int   Size();
+		void  clear();
+		void  Configure(Elements* elements,Loads* loads, DataSet* nodes, Vertices* vertices, Materials* materials,Parameters* parameters);
+		Object* GetObjectByOffset(int offset);
+		Object* GetObjectById(int* poffset,int eid);
+		void  Presort();
+		void  SetSorting(int* in_sorted_ids,int* in_id_offsets);
+		void  Sort();
+		DataSet* Copy(void);
+		int   DeleteObject(Object* object);
+		/*}}}*/
+
+};
+
+/*This routine cannot be object oriented, but need for demarshalling: */
+DataSet* DataSetDemarshall(char* marshalled_dataset);
+DataSet* DataSetDemarshallRaw(char** pmarshalled_dataset);
+
+
+/********************************************************INPUT************************************************/
+
+class Input;
+class Node;
+
+class Inputs: public DataSet{
+
+	public:
+
+		/*constructors, destructors: {{{1*/
+		Inputs();
+		Inputs(int enum_type);
+		~Inputs();
+		/*}}}*/
+		/*numerics: {{{1*/
+		int  AddInput(Input* in_input);
+		Input* GetInput(int enum_name);
+		Inputs* SpawnTriaInputs(int* indices);
+		Inputs* SpawnBeamInputs(int* indices);
+		Inputs* SpawnSingInputs(int  index  );
+		
+		void GetParameterValue(bool* pvalue,int enum_type);
+		void GetParameterValue(int* pvalue,int enum_type);
+		void GetParameterValue(double* pvalue,int enum_type);
+		void GetParameterValue(double* pvalue,Node* node,int enum_type);
+		void GetParameterValue(double* pvalue,Node* node1,Node* node2,double gauss_coord,int enum_type);
+		void GetParameterValue(double* pvalue,double* gauss,int enum_type);
+		void GetParameterValue(double* pvalue,double* gauss,int enum_type,double defaultvalue);
+		void GetParameterAverage(double* pvalue, int enum_type);
+		
+		void GetParameterValues(double* values,double* gauss_pointers, int numgauss,int enum_type);
+		void GetParameterValues(double* values,double* gauss_pointers, int numgauss,int enum_type,double* defaultvalues);
+	
+		void GetParameterDerivativeValue(double* derivativevalues, double* xyz_list, double* gauss,int enum_type);
+		void GetStrainRate2d(double* epsilon,double* xyz_list, double* gauss, int vxenum, int vyenum);
+		void GetStrainRate3d(double* epsilon,double* xyz_list, double* gauss, int vxenum, int vyenum,int vzenum);
+		void GetStrainRate3dPattyn(double* epsilon,double* xyz_list, double* gauss, int vxenum, int vyenum);
+
+		void ChangeEnum(int enumtype,int new_enumtype);
+
+		/*}}}*/
+
+};
+
+
+/********************************************************NODES************************************************/
+
+class Nodes: public DataSet{
+
+	public:
+
+		/*constructors, destructors: {{{1*/
+		Nodes();
+		Nodes(int enum_type);
+		~Nodes();
+		/*}}}*/
+		/*numerics: {{{1*/
+		void  DistributeDofs(int numberofnodes,int numdofspernode,int analysis_type);
+		void  FlagClones(int numberofnodes,int analysis_type);
+		void  FlagNodeSets(Vec pv_g, Vec pv_m, Vec pv_n, Vec pv_f, Vec pv_s,int analysis_type);
+		int   NumberOfDofs(int analysis_type);
+		int   NumberOfNodes(int analysis_type);
+		int   NumberOfNodes(void);
+		void  Ranks(int* ranks,int analysis_type);
+		/*}}}*/
+
+};
+
+
+/********************************************************VERTICES************************************************/
+
+class Vertices: public DataSet{
+
+	public:
+
+		/*constructors, destructors: {{{1*/
+		Vertices();
+		Vertices(int enum_type);
+		~Vertices();
+		/*}}}*/
+		/*numerics: {{{1*/
+		void  CreatePartitioningVector(Vec* ppartition,int numobjects);
+		void  DistributeDofs(int numberofnodes,int numdofspernode);
+		void  FlagClones(int numberofnodes);
+		int   NumberOfVertices(void);
+		void  Ranks(int* ranks);
+		/*}}}*/
+
+};
+
+
+/********************************************************LOADS************************************************/
+
+class Loads: public DataSet{
+
+	public:
+
+		/*constructors, destructors: {{{1*/
+		Loads();
+		Loads(int enum_type);
+		~Loads();
+		/*}}}*/
+		/*numerics: {{{1*/
+		int   MeltingIsPresent();
+		void  MeltingConstraints(int* pconverged, int* pnum_unstable_constraints);
+		int   NumberOfLoads(void);
+		void  OutputRifts(Vec riftproperties);
+		/*}}}*/
+
+};
+
+/********************************************************ELEMENTS************************************************/
+
+class Elements: public DataSet{
+
+	public:
+
+		/*constructors, destructors: {{{1*/
+		Elements();
+		Elements(int enum_type);
+		~Elements();
+		/*}}}*/
+		/*numerics: {{{1*/
+		/*}}}*/
+
+};
+
+/********************************************************CONSTRAINTS************************************************/
+
+class Constraints: public DataSet{
+
+	public:
+
+		/*constructors, destructors: {{{1*/
+		Constraints();
+		Constraints(int enum_type);
+		~Constraints();
+		/*}}}*/
+		/*numerics: {{{1*/
+		int   NumberOfLocalRgbs(int analysis_type);
+		int   NumberOfConstraints(void);
+		void  SetupSpcs(Nodes* nodes,Vec yg,int analysis_type);
+		void  SetupMpcs(Mat Rmg,Nodes* nodes,int analysis_type);
+		/*}}}*/
+
+};
+
+/********************************************************MATERIALS************************************************/
+
+class Materials: public DataSet{
+
+	public:
+
+		/*constructors, destructors: {{{1*/
+		Materials();
+		Materials(int enum_type);
+		~Materials();
+		/*}}}*/
+		/*numerics: {{{1*/
+		/*}}}*/
+
+};
+
+
+/********************************************************PARAMETERS************************************************/
+
+class Parameters: public DataSet{
+
+	public:
+
+		/*constructors, destructors: {{{1*/
+		Parameters();
+		Parameters(int enum_type);
+		~Parameters();
+		/*}}}*/
+		/*numerics: {{{1*/
+		int   FindParam(bool* pinteger,int enum_type);
+		int   FindParam(int* pinteger,int enum_type);
+		int   FindParam(double* pscalar, int enum_type);
+		int   FindParam(char** pstring,int enum_type);
+		int   FindParam(char*** pstringarray,int* pM,int enum_type);
+		int   FindParam(double** pdoublearray,int* pM,int enum_type);
+		int   FindParam(double** pdoublearray,int* pM,int* pN,int enum_type);
+		int   FindParam(Vec* pvec,int enum_type);
+		int   FindParam(Mat* pmat,int enum_type);
+		
+		void  SetParam(bool boolean,int enum_type);
+		void  SetParam(int integer,int enum_type);
+		void  SetParam(double scalar, int enum_type);
+		void  SetParam(char* string,int enum_type);
+		void  SetParam(char** stringarray,int M,int enum_type);
+		void  SetParam(double* doublearray,int M,int enum_type);
+		void  SetParam(double* doublearray,int M,int N,int enum_type);
+		void  SetParam(Vec vec,int enum_type);
+		void  SetParam(Mat mat,int enum_type);
+
+		Object* FindParamObject(int enum_type);
+		/*}}}*/
+
+};
+
+
+/********************************************************RESULTS************************************************/
+
+class Results: public DataSet{
+
+	public:
+
+		/*constructors, destructors: {{{1*/
+		Results();
+		Results(int enum_type);
+		~Results();
+		/*}}}*/
+		/*numerics: {{{1*/
+		Results* SpawnTriaResults(int* indices);
+		Results* SpawnBeamResults(int* indices);
+		Results* SpawnSingResults(int  index  );
+		/*}}}*/
+
+};
+
+
+
+#endif
Index: /issm/trunk/src/c/Container/Elements.cpp
===================================================================
--- /issm/trunk/src/c/Container/Elements.cpp	(revision 4235)
+++ /issm/trunk/src/c/Container/Elements.cpp	(revision 4235)
@@ -0,0 +1,44 @@
+/*
+ * \file Elements.c
+ * \brief: implementation of the Elements class, derived from DataSet class
+ */
+
+/*Headers: {{{1*/
+#ifdef HAVE_CONFIG_H
+	#include "config.h"
+#else
+#error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!"
+#endif
+
+#include <vector>
+#include <functional>
+#include <algorithm>
+#include <iostream>
+
+#include "./DataSet.h"
+#include "../shared/shared.h"
+#include "../include/include.h"
+#include "../EnumDefinitions/EnumDefinitions.h"
+
+using namespace std;
+/*}}}*/
+
+/*Object constructors and destructor*/
+/*FUNCTION Elements::Elements(){{{1*/
+Elements::Elements(){
+	return;
+}
+/*}}}*/
+/*FUNCTION Elements::Elements(int in_enum){{{1*/
+Elements::Elements(int in_enum): DataSet(in_enum){
+	//do nothing;
+	return;
+}
+/*}}}*/
+/*FUNCTION Elements::~Elements(){{{1*/
+Elements::~Elements(){
+	return;
+}
+/*}}}*/
+
+/*Object management*/
Index: /issm/trunk/src/c/Container/Inputs.cpp
===================================================================
--- /issm/trunk/src/c/Container/Inputs.cpp	(revision 4235)
+++ /issm/trunk/src/c/Container/Inputs.cpp	(revision 4235)
@@ -0,0 +1,655 @@
+/*
+ * \file Inputs.c
+ * \brief: implementation of the Inputs class, derived from DataSet class
+ */
+
+/*Headers: {{{1*/
+#ifdef HAVE_CONFIG_H
+	#include "config.h"
+#else
+#error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!"
+#endif
+
+#include <vector>
+#include <functional>
+#include <algorithm>
+#include <iostream>
+
+#include "./DataSet.h"
+#include "../shared/shared.h"
+#include "../include/include.h"
+#include "../EnumDefinitions/EnumDefinitions.h"
+
+using namespace std;
+/*}}}*/
+
+/*Object constructors and destructor*/
+/*FUNCTION Inputs::Inputs(){{{1*/
+Inputs::Inputs(){
+	return;
+}
+/*}}}*/
+/*FUNCTION Inputs::Inputs(int in_enum){{{1*/
+Inputs::Inputs(int in_enum): DataSet(in_enum) {
+	//do nothing;
+	return;
+}
+/*}}}*/
+/*FUNCTION Inputs::~Inputs(){{{1*/
+Inputs::~Inputs(){
+	return;
+}
+/*}}}*/
+
+/*Object management*/
+/*FUNCTION Inputs::GetParameterValue(double* pvalue,double* gauss,int enum_type){{{1*/
+void Inputs::GetParameterValue(double* pvalue,double* gauss, int enum_type){
+
+	vector<Object*>::iterator object;
+	Input* input=NULL;
+	bool   found=false;
+
+	/*Go through inputs and check whether any input with the same name is already in: */
+	for ( object=objects.begin() ; object < objects.end(); object++ ){
+
+		input=(Input*)(*object); 
+		if (input->EnumType()==enum_type){
+			found=true;
+			break;
+		}
+	}
+
+	if (!found){
+		/*we could not find an input with the correct enum type. No defaults values were provided, 
+		 * error out: */
+		ISSMERROR("could not find input with enum type %i (%s)",enum_type,EnumAsString(enum_type));
+	}
+
+	/*Ok, we have an input if we made it here, request the input to return the values: */
+	input->GetParameterValue(pvalue,gauss);
+
+}
+/*}}}*/
+/*FUNCTION Inputs::GetParameterValue(double* pvalue,double* gauss,int enum_type,double defaultvalue){{{1*/
+void Inputs::GetParameterValue(double* pvalue,double* gauss, int enum_type,double defaultvalue){
+
+	vector<Object*>::iterator object;
+	Input* input=NULL;
+	bool   found=false;
+
+	/*Go through inputs and check whether any input with the same name is already in: */
+	for ( object=objects.begin() ; object < objects.end(); object++ ){
+
+		input=(Input*)(*object); 
+		if (input->EnumType()==enum_type){
+			found=true;
+			break;
+		}
+	}
+
+	if (!found){
+		/*we could not find an input with the correct enum type. Return the default value: */
+		*pvalue=defaultvalue;
+	}
+	else{
+		input->GetParameterValue(pvalue,gauss);
+	}
+}
+/*}}}*/
+/*FUNCTION Inputs::GetParameterValues(double* values,double* gauss_pointers, int numgauss,int enum_type){{{1*/
+void Inputs::GetParameterValues(double* values,double* gauss_pointers, int numgauss,int enum_type){
+
+	vector<Object*>::iterator object;
+	Input* input=NULL;
+	bool   found=false;
+
+	/*Go through inputs and check whether any input with the same name is already in: */
+	for ( object=objects.begin() ; object < objects.end(); object++ ){
+
+		input=(Input*)(*object); 
+		if (input->EnumType()==enum_type){
+			found=true;
+			break;
+		}
+	}
+
+	if (!found){
+		/*we could not find an input with the correct enum type. No defaults values were provided, 
+		 * error out: */
+		ISSMERROR("could not find input with enum type %i (%s)",enum_type,EnumAsString(enum_type));
+	}
+
+	/*Ok, we have an input if we made it here, request the input to return the values: */
+	input->GetParameterValues(values,gauss_pointers,numgauss);
+
+}
+/*}}}*/
+/*FUNCTION Inputs::GetParameterValue(double* pvalue, Node* node, int enum_type){{{1*/
+void Inputs::GetParameterValue(double* pvalue,Node* node,int enum_type){
+
+	/*given a node, instead of a gauss point, we want to recover a value: probably in an element!: */
+
+	vector<Object*>::iterator object;
+	Input* input=NULL;
+	bool   found=false;
+
+	/*Go through inputs and check whether any input with the same name is already in: */
+	for ( object=objects.begin() ; object < objects.end(); object++ ){
+
+		input=(Input*)(*object); 
+		if (input->EnumType()==enum_type){
+			found=true;
+			break;
+		}
+	}
+
+	if (!found){
+		/*we could not find an input with the correct enum type. No defaults values were provided, 
+		 * error out: */
+		ISSMERROR("could not find input with enum type %i (%s)",enum_type,EnumAsString(enum_type));
+	}
+
+	/*Ok, we have an input if we made it here, request the input to return the values: */
+	input->GetParameterValue(pvalue,node);
+}
+/*}}}*/
+/*FUNCTION Inputs::GetParameterValue(double* pvalue, Node* node1, Node* node2,int enum_type){{{1*/
+void Inputs::GetParameterValue(double* pvalue,Node* node1, Node* node2,double gauss_coord,int enum_type){
+
+	/*given a node, instead of a gauss point, we want to recover a value: probably in an element!: */
+
+	vector<Object*>::iterator object;
+	Input* input=NULL;
+	bool   found=false;
+
+	/*Go through inputs and check whether any input with the same name is already in: */
+	for ( object=objects.begin() ; object < objects.end(); object++ ){
+
+		input=(Input*)(*object); 
+		if (input->EnumType()==enum_type){
+			found=true;
+			break;
+		}
+	}
+
+	if (!found){
+		/*we could not find an input with the correct enum type. No defaults values were provided, 
+		 * error out: */
+		ISSMERROR("could not find input with enum type %i (%s)",enum_type,EnumAsString(enum_type));
+	}
+
+	/*Ok, we have an input if we made it here, request the input to return the values: */
+	input->GetParameterValue(pvalue,node1,node2,gauss_coord);
+}
+/*}}}*/
+/*FUNCTION Inputs::GetParameterValues(double* values,double* gauss_pointers, int numgauss,int enum_type,double* defaultvalues){{{1*/
+void Inputs::GetParameterValues(double* values,double* gauss_pointers, int numgauss,int enum_type,double* defaultvalues){
+
+	int i;
+	vector<Object*>::iterator object;
+	Input* input=NULL;
+	bool   found=false;
+
+	/*Go through inputs and check whether any input with the same name is already in: */
+	for ( object=objects.begin() ; object < objects.end(); object++ ){
+
+		input=(Input*)(*object); 
+		if (input->EnumType()==enum_type){
+			found=true;
+			break;
+		}
+	}
+
+	if (!found){
+		/*we could not find an input with the correct enum type. Return the default values: */
+		for(i=0;i<numgauss;i++) values[i]=defaultvalues[i];
+	}
+	else{
+		input->GetParameterValues(values,gauss_pointers,numgauss);
+	}
+
+}
+/*}}}*/
+/*FUNCTION Inputs::GetParameterAverage(double* pvalue,int enum-type){{{1*/
+void Inputs::GetParameterAverage(double* pvalue,int enum_type){
+
+	vector<Object*>::iterator object;
+	Input* input=NULL;
+	bool   found=false;
+
+	/*Go through inputs and check whether any input with the same name is already in: */
+	for ( object=objects.begin() ; object < objects.end(); object++ ){
+
+		input=(Input*)(*object); 
+		if (input->EnumType()==enum_type){
+			found=true;
+			break;
+		}
+	}
+
+	if (!found){
+		/*we could not find an input with the correct enum type. No defaults values were provided, 
+		 * error out: */
+		ISSMERROR("could not find input with enum type %i (%s)",enum_type,EnumAsString(enum_type));
+	}
+
+	/*Ok, we have an input if we made it here, request the input to return the value: */
+	input->GetParameterAverage(pvalue);
+
+}
+/*}}}*/
+/*FUNCTION Inputs::GetParameterValue(bool* pvalue,int enum-type){{{1*/
+void Inputs::GetParameterValue(bool* pvalue,int enum_type){
+
+	vector<Object*>::iterator object;
+	Input* input=NULL;
+	bool   found=false;
+
+	/*Go through inputs and check whether any input with the same name is already in: */
+	for ( object=objects.begin() ; object < objects.end(); object++ ){
+
+		input=(Input*)(*object); 
+		if (input->EnumType()==enum_type){
+			found=true;
+			break;
+		}
+	}
+
+	if (!found){
+		/*we could not find an input with the correct enum type. No defaults values were provided, 
+		 * error out: */
+		ISSMERROR("could not find input with enum type %i (%s)",enum_type,EnumAsString(enum_type));
+	}
+
+	/*Ok, we have an input if we made it here, request the input to return the value: */
+	input->GetParameterValue(pvalue);
+
+}
+/*}}}*/
+/*FUNCTION Inputs::GetParameterValue(int* pvalue,int enum-type){{{1*/
+void Inputs::GetParameterValue(int* pvalue,int enum_type){
+
+	vector<Object*>::iterator object;
+	Input* input=NULL;
+	bool   found=false;
+
+	/*Go through inputs and check whether any input with the same name is already in: */
+	for ( object=objects.begin() ; object < objects.end(); object++ ){
+
+		input=(Input*)(*object); 
+		if (input->EnumType()==enum_type){
+			found=true;
+			break;
+		}
+	}
+
+	if (!found){
+		/*we could not find an input with the correct enum type. No defaults values were provided, 
+		 * error out: */
+		ISSMERROR("could not find input with enum type %i (%s)",enum_type,EnumAsString(enum_type));
+	}
+
+	/*Ok, we have an input if we made it here, request the input to return the value: */
+	input->GetParameterValue(pvalue);
+
+}
+/*}}}*/
+/*FUNCTION Inputs::GetParameterValue(double* pvalue,int enum-type){{{1*/
+void Inputs::GetParameterValue(double* pvalue,int enum_type){
+
+	vector<Object*>::iterator object;
+	Input* input=NULL;
+	bool   found=false;
+
+	/*Go through inputs and check whether any input with the same name is already in: */
+	for ( object=objects.begin() ; object < objects.end(); object++ ){
+
+		input=(Input*)(*object); 
+		if (input->EnumType()==enum_type){
+			found=true;
+			break;
+		}
+	}
+
+	if (!found){
+		/*we could not find an input with the correct enum type. No defaults values were provided, 
+		 * error out: */
+		ISSMERROR("could not find input with enum type %i (%s)",enum_type,EnumAsString(enum_type));
+	}
+
+	/*Ok, we have an input if we made it here, request the input to return the value: */
+	input->GetParameterValue(pvalue);
+
+}
+/*}}}*/
+/*FUNCTION Inputs::GetParameterDerivativeValue(double* derivativevalues, double* xyz_list, double* gauss,int enum_type){{{1*/
+void Inputs::GetParameterDerivativeValue(double* derivativevalues, double* xyz_list, double* gauss,int enum_type){
+
+	vector<Object*>::iterator object;
+	Input* input=NULL;
+	bool   found=false;
+
+	/*Go through inputs and check whether any input with the same name is already in: */
+	for ( object=objects.begin() ; object < objects.end(); object++ ){
+
+		input=(Input*)(*object); 
+		if (input->EnumType()==enum_type){
+			found=true;
+			break;
+		}
+	}
+
+	if (!found){
+		/*we could not find an input with the correct enum type. No defaults values were provided, 
+		 * error out: */
+		ISSMERROR("could not find input with enum type %i (%s)",enum_type,EnumAsString(enum_type));
+	}
+
+	/*Ok, we have an input if we made it here, request the input to return the value: */
+	input->GetParameterDerivativeValue(derivativevalues,xyz_list,gauss);
+}
+/*}}}*/
+/*FUNCTION Inputs::GetStrainRate2d(double* epsilon,double* xyz_list, double* gauss, int vxenum, int vyenum){{{1*/
+void Inputs::GetStrainRate2d(double* epsilon,double* xyz_list, double* gauss, int vxenum, int vyenum){
+	/*Compute the 2d Strain Rate (3 components):
+	 *
+	 * epsilon=[exx eyy exy]
+	 */
+
+	vector<Object*>::iterator object;
+	int i;
+	Input* vxinput=NULL;
+	Input* vyinput=NULL;
+	double epsilonvx[3];
+	double epsilonvy[3];
+	bool   foundvx=false;
+	bool   foundvy=false;
+
+	/*Go through inputs and find data for vxenum: */
+	for ( object=objects.begin() ; object < objects.end(); object++ ){
+		vxinput=(Input*)(*object); 
+		if (vxinput->EnumType()==vxenum){
+			foundvx=true;
+			break;
+		}
+	}
+	/*Go through inputs and find data for vyenum: */
+	for ( object=objects.begin() ; object < objects.end(); object++ ){
+		vyinput=(Input*)(*object); 
+		if (vyinput->EnumType()==vyenum){
+			foundvy=true;
+			break;
+		}
+	}
+
+	/*Check that both inputs have been found*/
+	if (!foundvx || !foundvy){
+		ISSMERROR("Could not find input with enum %i (%s) or enum %i (%s)",vxenum,EnumAsString(vxenum),vyenum,EnumAsString(vyenum));
+	}
+
+	/*Get strain rate assuming that epsilon has been allocated*/
+	vxinput->GetVxStrainRate2d(epsilonvx,xyz_list,gauss);
+	vyinput->GetVyStrainRate2d(epsilonvy,xyz_list,gauss);
+
+	/*Sum all contributions*/
+	for(i=0;i<3;i++) epsilon[i]=epsilonvx[i]+epsilonvy[i];
+
+}
+/*}}}*/
+/*FUNCTION Inputs::GetStrainRate3d(double* epsilon,double* xyz_list, double* gauss, int vxenum, int vyenum, int vzenum){{{1*/
+void Inputs::GetStrainRate3d(double* epsilon,double* xyz_list, double* gauss, int vxenum, int vyenum,int vzenum){
+	/*Compute the 3d Strain Rate (6 components):
+	 *
+	 * epsilon=[exx eyy ezz exy exz euz]
+	 */
+
+	int    i;
+	vector<Object*>::iterator object;
+	Input* vxinput=NULL;
+	Input* vyinput=NULL;
+	Input* vzinput=NULL;
+	bool   foundvx=false;
+	bool   foundvy=false;
+	bool   foundvz=false;
+	double epsilonvx[6];
+	double epsilonvy[6];
+	double epsilonvz[6];
+
+	/*Go through inputs and find data for vxenum: */
+	for ( object=objects.begin() ; object < objects.end(); object++ ){
+		vxinput=(Input*)(*object); 
+		if (vxinput->EnumType()==vxenum){
+			foundvx=true;
+			break;
+		}
+	}
+	/*Go through inputs and find data for vyenum: */
+	for ( object=objects.begin() ; object < objects.end(); object++ ){
+		vyinput=(Input*)(*object); 
+		if (vyinput->EnumType()==vyenum){
+			foundvy=true;
+			break;
+		}
+	}
+	/*Go through inputs and find data for vzenum, not for Pattyn*/
+	for ( object=objects.begin() ; object < objects.end(); object++ ){
+		vzinput=(Input*)(*object); 
+		if (vzinput->EnumType()==vzenum){
+			foundvz=true;
+			break;
+		}
+	}
+
+	/*Check that all inputs have been found*/
+	if (!foundvx || !foundvy || !foundvz){
+		ISSMERROR("Could not find input with enum %i (%s), enum %i (%s) or  enum %i (%s)",vxenum,EnumAsString(vxenum),vyenum,EnumAsString(vyenum),vzenum,EnumAsString(vzenum));
+	}
+
+	/*Get strain rate assuming that epsilon has been allocated*/
+	vxinput->GetVxStrainRate3d(epsilonvx,xyz_list,gauss);
+	vyinput->GetVyStrainRate3d(epsilonvy,xyz_list,gauss);
+	vzinput->GetVzStrainRate3d(epsilonvz,xyz_list,gauss);
+
+	/*Sum all contributions*/
+	for(i=0;i<6;i++) epsilon[i]=epsilonvx[i]+epsilonvy[i]+epsilonvz[i];
+
+}
+/*}}}*/
+/*FUNCTION Inputs::GetStrainRate3dPattyn(double* epsilon,double* xyz_list, double* gauss, int vxenum, int vyenum){{{1*/
+void Inputs::GetStrainRate3dPattyn(double* epsilon,double* xyz_list, double* gauss, int vxenum, int vyenum){
+	/*Compute the 3d Blatter/PattynStrain Rate (5 components):
+	 *
+	 * epsilon=[exx eyy exy exz euz]
+	 *
+	 * with exz=1/2 du/dz
+	 *      eyz=1/2 dv/dz
+	 *
+	 * the contribution of vz is neglected
+	 */
+
+	int    i;
+	vector<Object*>::iterator object;
+	Input* vxinput=NULL;
+	Input* vyinput=NULL;
+	bool   foundvx=false;
+	bool   foundvy=false;
+	double epsilonvx[5];
+	double epsilonvy[5];
+
+	/*Go through inputs and find data for vxenum: */
+	for ( object=objects.begin() ; object < objects.end(); object++ ){
+		vxinput=(Input*)(*object); 
+		if (vxinput->EnumType()==vxenum){
+			foundvx=true;
+			break;
+		}
+	}
+	/*Go through inputs and find data for vyenum: */
+	for ( object=objects.begin() ; object < objects.end(); object++ ){
+		vyinput=(Input*)(*object); 
+		if (vyinput->EnumType()==vyenum){
+			foundvy=true;
+			break;
+		}
+	}
+
+	/*Check that all inputs have been found*/
+	if (!foundvx || !foundvy){
+		ISSMERROR("Could not find input with enum %i (%s) or enum %i (%s)",vxenum,EnumAsString(vxenum),vyenum,EnumAsString(vyenum));
+	}
+
+	/*Get strain rate assuming that epsilon has been allocated*/
+	vxinput->GetVxStrainRate3dPattyn(epsilonvx,xyz_list,gauss);
+	vyinput->GetVyStrainRate3dPattyn(epsilonvy,xyz_list,gauss);
+
+	/*Sum all contributions*/
+	for(i=0;i<5;i++) epsilon[i]=epsilonvx[i]+epsilonvy[i];
+
+}
+/*}}}*/
+/*FUNCTION Inputs::AddInput{{{1*/
+int  Inputs::AddInput(Input* in_input){
+
+	/*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: */
+	vector<Object*>::iterator object;
+	Input* input=NULL;
+
+	for ( object=objects.begin() ; object < objects.end(); object++ ){
+
+		input=(Input*)(*object); 
+
+		if (input->EnumType()==in_input->EnumType()){
+			this->DeleteObject(input);
+			break;
+		}
+	}
+	this->AddObject(in_input);
+}
+/*}}}*/
+/*FUNCTION Inputs::GetInput{{{1*/
+Input* Inputs::GetInput(int enum_name){
+
+	vector<Object*>::iterator object;
+	Input* input=NULL;
+
+	for ( object=objects.begin() ; object < objects.end(); object++ ){
+
+		input=(Input*)(*object); 
+
+		if (input->EnumType()==enum_name){
+			return input;
+		}
+	}
+	return NULL;
+}
+/*}}}*/
+/*FUNCTION Inputs::ChangeEnum{{{1*/
+void  Inputs::ChangeEnum(int oldenumtype,int newenumtype){
+
+	/*Go through dataset of inputs and look for input with 
+	 * same enum as input enum, once found, just change its name */
+	vector<Object*>::iterator object;
+	Input* input=NULL;
+
+	/*Delete existing input of newenumtype if it exists*/
+	for ( object=objects.begin() ; object < objects.end(); object++ ){
+		input=(Input*)(*object); 
+
+		if (input->EnumType()==newenumtype){
+			this->DeleteObject(input);
+			break;
+		}
+	}
+
+	/*Change enum_type of input of oldenumtype*/
+	for ( object=objects.begin() ; object < objects.end(); object++ ){
+
+		input=(Input*)(*object); 
+
+		if (input->EnumType()==oldenumtype){
+			input->ChangeEnum(newenumtype);
+			break;
+		}
+	}
+}
+/*}}}*/
+/*FUNCTION Inputs::SpawnBeamInputs{{{1*/
+Inputs* Inputs::SpawnBeamInputs(int* indices){
+
+	/*Intermediary*/
+	vector<Object*>::iterator object;
+	Input* inputin=NULL;
+	Input* inputout=NULL;
+
+	/*Output*/
+	Inputs* newinputs=new Inputs();
+
+	/*Go through inputs and call Spawn function*/
+	for ( object=objects.begin() ; object < objects.end(); object++ ){
+
+		/*Create new input*/
+		inputin=(Input*)(*object); 
+		inputout=inputin->SpawnBeamInput(indices);
+
+		/*Add input to new inputs*/
+		newinputs->AddObject(inputout);
+	}
+
+	/*Assign output pointer*/
+	return newinputs;
+}
+/*}}}*/
+/*FUNCTION Inputs::SpawnSingInputs{{{1*/
+Inputs* Inputs::SpawnSingInputs(int index){
+
+	/*Intermediary*/
+	vector<Object*>::iterator object;
+	Input* inputin=NULL;
+	Input* inputout=NULL;
+
+	/*Output*/
+	Inputs* newinputs=new Inputs();
+
+	/*Go through inputs and call Spawn function*/
+	for ( object=objects.begin() ; object < objects.end(); object++ ){
+
+		/*Create new input*/
+		inputin=(Input*)(*object); 
+		inputout=inputin->SpawnSingInput(index);
+
+		/*Add input to new inputs*/
+		newinputs->AddObject(inputout);
+	}
+
+	/*Assign output pointer*/
+	return newinputs;
+}
+/*}}}*/
+/*FUNCTION Inputs::SpawnTriaInputs{{{1*/
+Inputs* Inputs::SpawnTriaInputs(int* indices){
+
+	/*Intermediary*/
+	vector<Object*>::iterator object;
+	Input* inputin=NULL;
+	Input* inputout=NULL;
+
+	/*Output*/
+	Inputs* newinputs=new Inputs();
+
+	/*Go through inputs and call Spawn function*/
+	for ( object=objects.begin() ; object < objects.end(); object++ ){
+
+		/*Create new input*/
+		inputin=(Input*)(*object); 
+		inputout=inputin->SpawnTriaInput(indices);
+
+		/*Add input to new inputs*/
+		newinputs->AddObject(inputout);
+	}
+
+	/*Assign output pointer*/
+	return newinputs;
+}
+/*}}}*/
Index: /issm/trunk/src/c/Container/Loads.cpp
===================================================================
--- /issm/trunk/src/c/Container/Loads.cpp	(revision 4235)
+++ /issm/trunk/src/c/Container/Loads.cpp	(revision 4235)
@@ -0,0 +1,143 @@
+/*
+ * \file Loads.c
+ * \brief: implementation of the Loads class, derived from DataSet class
+ */
+
+/*Headers: {{{1*/
+#ifdef HAVE_CONFIG_H
+	#include "config.h"
+#else
+#error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!"
+#endif
+
+#include <vector>
+#include <functional>
+#include <algorithm>
+#include <iostream>
+
+#include "./DataSet.h"
+#include "../shared/shared.h"
+#include "../include/include.h"
+#include "../EnumDefinitions/EnumDefinitions.h"
+
+using namespace std;
+/*}}}*/
+
+/*Object constructors and destructor*/
+/*FUNCTION Loads::Loads(){{{1*/
+Loads::Loads(){
+	return;
+}
+/*}}}*/
+/*FUNCTION Loads::Loads(int in_enum){{{1*/
+Loads::Loads(int in_enum): DataSet(in_enum){
+	//do nothing;
+	return;
+}
+/*}}}*/
+/*FUNCTION Loads::~Loads(){{{1*/
+Loads::~Loads(){
+	return;
+}
+/*}}}*/
+
+/*Numerics:*/
+/*FUNCTION Loads::MeltingIsPresent{{{1*/
+int   Loads::MeltingIsPresent(){
+
+	int i;
+	int found=0;
+	int mpi_found=0;
+
+	for(i=0;i<this->Size();i++){
+		Object* object=(Object*)this->GetObjectByOffset(i);
+		if (object->Enum()==PengridEnum){
+			found=1;
+			break;
+		}
+	}
+	
+	#ifdef _PARALLEL_
+	MPI_Reduce (&found,&mpi_found,1,MPI_INT,MPI_SUM,0,MPI_COMM_WORLD );
+	MPI_Bcast(&mpi_found,1,MPI_INT,0,MPI_COMM_WORLD);                
+	found=mpi_found;
+	#endif
+
+	return found;
+}
+/*}}}*/
+/*FUNCTION Loads::MeltingConstraints{{{1*/
+void  Loads::MeltingConstraints(int* pconverged, int* pnum_unstable_constraints){
+
+	int i;
+
+	int unstable=0;
+	int num_unstable_constraints=0;
+	int converged=0;
+	int sum_num_unstable_constraints=0;
+
+	num_unstable_constraints=0;	
+
+	/*Enforce constraints: */
+	for(i=0;i<this->Size();i++){
+		Object* object=(Object*)this->GetObjectByOffset(i);
+		if(object->Enum()==PengridEnum){
+
+			Pengrid* pengrid=(Pengrid*)object;
+			pengrid->PenaltyConstrain(&unstable);
+			num_unstable_constraints+=unstable;
+		}
+	}
+
+	#ifdef _PARALLEL_
+	MPI_Reduce (&num_unstable_constraints,&sum_num_unstable_constraints,1,MPI_INT,MPI_SUM,0,MPI_COMM_WORLD );
+	MPI_Bcast(&sum_num_unstable_constraints,1,MPI_INT,0,MPI_COMM_WORLD);                
+	num_unstable_constraints=sum_num_unstable_constraints;
+	#endif
+
+	/*Have we converged? : */
+	if (num_unstable_constraints==0) converged=1;
+	else converged=0;
+
+	/*Assign output pointers: */
+	*pconverged=converged;
+	*pnum_unstable_constraints=num_unstable_constraints;
+}
+/*}}}*/
+/*FUNCTION Loads::NumberOfLoads{{{1*/
+int Loads::NumberOfLoads(void){
+
+	int localloads;
+	int numberofloads;
+
+	/*Get number of local loads*/
+	localloads=this->Size();
+
+	/*figure out total number of loads combining all the cpus (no clones here)*/
+	#ifdef _PARALLEL_
+	MPI_Reduce(&localloads,&numberofloads,1,MPI_INT,MPI_SUM,0,MPI_COMM_WORLD );
+	MPI_Bcast(&numberofloads,1,MPI_INT,0,MPI_COMM_WORLD);
+	#else
+	numberofloads=localloads;
+	#endif
+
+	return numberofloads;
+}
+/*}}}*/
+/*FUNCTION Loads::OutputRifts{{{1*/
+void  Loads::OutputRifts(Vec riftproperties){
+
+	int i;
+
+	for(i=0;i<this->Size();i++){
+
+		Object* object=(Object*)this->GetObjectByOffset(i);
+
+		if(object->Enum()==RiftfrontEnum){
+
+			Riftfront* riftfront=(Riftfront*)object;
+			riftfront->OutputProperties(riftproperties);
+		}
+	}
+}
+/*}}}*/
Index: /issm/trunk/src/c/Container/Materials.cpp
===================================================================
--- /issm/trunk/src/c/Container/Materials.cpp	(revision 4235)
+++ /issm/trunk/src/c/Container/Materials.cpp	(revision 4235)
@@ -0,0 +1,44 @@
+/*
+ * \file Materials.c
+ * \brief: implementation of the Materials class, derived from DataSet class
+ */
+
+/*Headers: {{{1*/
+#ifdef HAVE_CONFIG_H
+	#include "config.h"
+#else
+#error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!"
+#endif
+
+#include <vector>
+#include <functional>
+#include <algorithm>
+#include <iostream>
+
+#include "./DataSet.h"
+#include "../shared/shared.h"
+#include "../include/include.h"
+#include "../EnumDefinitions/EnumDefinitions.h"
+
+using namespace std;
+/*}}}*/
+
+/*Object constructors and destructor*/
+/*FUNCTION Materials::Materials(){{{1*/
+Materials::Materials(){
+	return;
+}
+/*}}}*/
+/*FUNCTION Materials::Materials(int in_enum){{{1*/
+Materials::Materials(int in_enum): DataSet(in_enum){
+	//do nothing;
+	return;
+}
+/*}}}*/
+/*FUNCTION Materials::~Materials(){{{1*/
+Materials::~Materials(){
+	return;
+}
+/*}}}*/
+
+/*Object management*/
Index: /issm/trunk/src/c/Container/Nodes.cpp
===================================================================
--- /issm/trunk/src/c/Container/Nodes.cpp	(revision 4235)
+++ /issm/trunk/src/c/Container/Nodes.cpp	(revision 4235)
@@ -0,0 +1,329 @@
+/*
+ * \file Nodes.c
+ * \brief: implementation of the Nodes class, derived from DataSet class
+ */
+
+/*Headers: {{{1*/
+#ifdef HAVE_CONFIG_H
+	#include "config.h"
+#else
+#error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!"
+#endif
+
+#include <vector>
+#include <functional>
+#include <algorithm>
+#include <iostream>
+
+#include "./DataSet.h"
+#include "../shared/shared.h"
+#include "../include/include.h"
+#include "../EnumDefinitions/EnumDefinitions.h"
+
+using namespace std;
+/*}}}*/
+
+/*Object constructors and destructor*/
+/*FUNCTION Nodes::Nodes(){{{1*/
+Nodes::Nodes(){
+	return;
+}
+/*}}}*/
+/*FUNCTION Nodes::Nodes(int in_enum){{{1*/
+Nodes::Nodes(int in_enum): DataSet(in_enum){
+	//do nothing;
+	return;
+}
+/*}}}*/
+/*FUNCTION Nodes::~Nodes(){{{1*/
+Nodes::~Nodes(){
+	return;
+}
+/*}}}*/
+
+/*Numerics*/
+/*FUNCTION Nodes::DistributeDofs{{{1*/
+void  Nodes::DistributeDofs(int numberofobjects,int numberofdofsperobject,int analysis_type){
+
+	extern int num_procs;
+	extern int my_rank;
+
+	int  i;
+	
+	int  dofcount=0;
+	int* alldofcount=NULL;
+	int* truedofs=NULL;
+	int* alltruedofs=NULL;
+
+	/*Go through objects, and distribute dofs locally, from 0 to numberofdofsperobject: */
+	for (i=0;i<this->Size();i++){
+		Node* node=(Node*)this->GetObjectByOffset(i);
+
+		/*Check that this node corresponds to our analysis currently being carried out: */
+		if (node->InAnalysis(analysis_type)){
+			node->DistributeDofs(&dofcount);
+		}
+	}
+
+	/*Ok, now every object has distributed dofs, but locally, and with a dof count starting from 
+	 *0. This means the dofs between all the cpus are not synchronized! We need to synchronize all 
+	 *dof on all cpus, and use those to update the dofs of every object: */
+
+	alldofcount=(int*)xmalloc(num_procs*sizeof(int));
+	MPI_Gather(&dofcount,1,MPI_INT,alldofcount,1,MPI_INT,0,MPI_COMM_WORLD);
+	MPI_Bcast(alldofcount,num_procs,MPI_INT,0,MPI_COMM_WORLD);
+
+	/*Ok, now every cpu should start its own dof count at the end of the dofcount 
+	 * from cpu-1. : */
+	dofcount=0;
+	if(my_rank==0){
+		dofcount=0;
+	}
+	else{
+		for(i=0;i<my_rank;i++){
+			dofcount+=alldofcount[i];
+		}
+	}
+
+
+	/*Ok, now every cpu knows where his dofs should start. Update the dof count: */
+	for (i=0;i<this->Size();i++){
+
+		/*Check that this node corresponds to our analysis currently being carried out: */
+		Node* node=(Node*)this->GetObjectByOffset(i);
+		if (node->InAnalysis(analysis_type)){
+			node->OffsetDofs(dofcount);
+		}
+
+	}
+
+	/*Finally, remember that cpus may have skipped some objects, because they were clones. For every 
+	 * object that is not a clone, tell them to show their dofs, so that later on, they can get picked 
+	 * up by their clones: */
+	truedofs=(int*)xcalloc(numberofobjects*numberofdofsperobject,sizeof(int));
+	alltruedofs=(int*)xcalloc(numberofobjects*numberofdofsperobject,sizeof(int));
+
+	for (i=0;i<this->Size();i++){
+		/*Check that this node corresponds to our analysis currently being carried out: */
+		Node* node=(Node*)this->GetObjectByOffset(i);
+		if (node->InAnalysis(analysis_type)){
+			node->ShowTrueDofs(truedofs);
+		}
+	}
+	
+	MPI_Allreduce ( (void*)truedofs,(void*)alltruedofs,numberofobjects*numberofdofsperobject,MPI_INT,MPI_MAX,MPI_COMM_WORLD);
+
+	/*Ok, now every cpu knows the true dofs of everyone else that is not a clone. Let the clones recover those true dofs: */
+	for (i=0;i<this->Size();i++){
+		/*Check that this node corresponds to our analysis currently being carried out: */
+		Node* node=(Node*)this->GetObjectByOffset(i);
+		if (node->InAnalysis(analysis_type)){
+			node->UpdateCloneDofs(alltruedofs);
+		}
+	}
+
+	/* Free ressources: */
+	xfree((void**)&alldofcount);
+	xfree((void**)&truedofs);
+	xfree((void**)&alltruedofs);
+
+}
+/*}}}*/
+/*FUNCTION Nodes::FlagClones{{{1*/
+void  Nodes::FlagClones(int numberofobjects,int analysis_type){
+
+	int i;
+	extern int num_procs;
+
+	int* ranks=NULL;
+	int* minranks=NULL;
+
+	/*Allocate ranks: */
+	ranks=(int*)xmalloc(numberofobjects*sizeof(int));
+	minranks=(int*)xmalloc(numberofobjects*sizeof(int));
+
+	for(i=0;i<numberofobjects;i++)ranks[i]=num_procs; //no cpu can have rank num_procs. This is the maximum limit.
+
+	/*Now go through all our objects and ask them to report to who they belong (which rank): */
+	Ranks(ranks,analysis_type);
+
+	/*We need to take the minimum rank for each vertex, and every cpu needs to get that result. That way, 
+	 * when we start building the dof list for all vertexs, a cpu can check whether its vertex already has been 
+	 * dealt with by another cpu. We take the minimum because we are going to manage dof assignment in increasing 
+	 * order of cpu rank. This is also why we initialized this array to num_procs.*/
+	MPI_Allreduce ( (void*)ranks,(void*)minranks,numberofobjects,MPI_INT,MPI_MIN,MPI_COMM_WORLD);
+
+	/*Now go through all objects, and use minranks to flag which objects are cloned: */
+	for(i=0;i<this->Size();i++){
+
+		Node* node=(Node*)this->GetObjectByOffset(i);
+
+		/*Check that this node corresponds to our analysis currently being carried out: */
+		if (node->InAnalysis(analysis_type)){
+
+			/*For this object, decide whether it is a clone: */
+			node->SetClone(minranks);
+		}
+	}
+
+	/*Free ressources: */
+	xfree((void**)&ranks); 
+	xfree((void**)&minranks);
+
+}
+/*}}}*/
+/*FUNCTION Nodes::FlagNodeSets{{{1*/
+void Nodes::FlagNodeSets(Vec pv_g, Vec pv_m, Vec pv_n, Vec pv_f, Vec pv_s,int analysis_type){
+
+	int i;
+
+	for(i=0;i<this->Size();i++){
+
+		Node* node=(Node*)this->GetObjectByOffset(i);
+			
+		/*Check that this node corresponds to our analysis currently being carried out: */
+		if (node->InAnalysis(analysis_type)){
+
+			/*Plug set values intp our 4 set vectors: */
+			node->CreateVecSets(pv_g,pv_m,pv_n,pv_f,pv_s);
+		}
+	}
+
+	/*Assemble: */
+	VecAssemblyBegin(pv_g);
+	VecAssemblyEnd(pv_g);
+
+	VecAssemblyBegin(pv_m);
+	VecAssemblyEnd(pv_m);
+
+	VecAssemblyBegin(pv_n);
+	VecAssemblyEnd(pv_n);
+
+	VecAssemblyBegin(pv_f);
+	VecAssemblyEnd(pv_f);
+
+	VecAssemblyBegin(pv_s);
+	VecAssemblyEnd(pv_s);
+
+}
+/*}}}*/
+/*FUNCTION Nodes::NumberOfDofs(int analysis_type){{{1*/
+int   Nodes::NumberOfDofs(int analysis_type){
+
+	int i;
+	
+	int   numdofs=0;
+	int   allnumdofs;
+
+	/*Now go through all nodes, and get how many dofs they own, unless they are clone nodes: */
+	for(i=0;i<this->Size();i++){
+
+		Node* node=(Node*)this->GetObjectByOffset(i);
+
+		/*Check that this node corresponds to our analysis currently being carried out: */
+		if (node->InAnalysis(analysis_type)){
+
+			/*Ok, this object is a node, ask it to plug values into partition: */
+			if (!node->IsClone()){
+
+				numdofs+=node->GetNumberOfDofs();
+
+			}
+		}
+	}
+
+	/*Gather from all cpus: */
+	MPI_Allreduce ( (void*)&numdofs,(void*)&allnumdofs,1,MPI_INT,MPI_SUM,MPI_COMM_WORLD);
+
+	return allnumdofs;
+}
+/*}}}*/
+/*FUNCTION Nodes::NumberOfNodes(){{{1*/
+int Nodes::NumberOfNodes(void){
+
+	int i;
+
+	int max_sid=0;
+	int sid;
+	int node_max_sid;
+
+	for(i=0;i<this->Size();i++){
+
+		Node* node=(Node*)this->GetObjectByOffset(i);
+		sid=node->Sid();
+		if (sid>max_sid)max_sid=sid;
+	}
+
+	#ifdef _PARALLEL_
+	MPI_Reduce (&max_sid,&node_max_sid,1,MPI_INT,MPI_MAX,0,MPI_COMM_WORLD );
+	MPI_Bcast(&node_max_sid,1,MPI_INT,0,MPI_COMM_WORLD);
+	max_sid=node_max_sid;
+	#endif 
+
+	/*sid starts at 0*/
+	max_sid++;
+
+	/*return*/
+	return max_sid;
+}
+/*}}}*/
+/*FUNCTION Nodes::NumberOfNodes{{{1*/
+int Nodes::NumberOfNodes(int analysis_type){
+
+	int i;
+
+	int max_sid=0;
+	int sid;
+	int node_max_sid;
+
+	for(i=0;i<this->Size();i++){
+
+		Node* node=(Node*)this->GetObjectByOffset(i);
+
+		/*Check that this node corresponds to our analysis currently being carried out: */
+		if (node->InAnalysis(analysis_type)){
+
+			sid=node->Sid();
+			if (sid>max_sid)max_sid=sid;
+		}
+	}
+
+#ifdef _PARALLEL_
+	MPI_Reduce (&max_sid,&node_max_sid,1,MPI_INT,MPI_MAX,0,MPI_COMM_WORLD );
+	MPI_Bcast(&node_max_sid,1,MPI_INT,0,MPI_COMM_WORLD);
+	max_sid=node_max_sid;
+#endif 
+
+	/*sid starts at 0*/
+	max_sid++;
+
+	/*return*/
+	return max_sid;
+}
+/*}}}*/
+/*FUNCTION Nodes::Ranks{{{1*/
+void   Nodes::Ranks(int* ranks,int analysis_type){
+
+	/*Go through nodes, and for each object, report it cpu: */
+
+	int i;
+	int rank;
+	int sid;
+	
+	for(i=0;i<this->Size();i++){
+
+		Node* node=(Node*)this->GetObjectByOffset(i);
+
+		/*Check that this node corresponds to our analysis currently being carried out: */
+		if (node->InAnalysis(analysis_type)){
+
+			rank=node->MyRank();
+			sid=node->Sid();
+
+			/*Plug rank into ranks, according to id: */
+			ranks[sid]=rank; 
+		}
+	}
+	return;
+}
+/*}}}*/
Index: /issm/trunk/src/c/Container/Parameters.cpp
===================================================================
--- /issm/trunk/src/c/Container/Parameters.cpp	(revision 4235)
+++ /issm/trunk/src/c/Container/Parameters.cpp	(revision 4235)
@@ -0,0 +1,416 @@
+/*
+ * \file Parameters.c
+ * \brief: implementation of the Parameters class, derived from DataSet class
+ */
+
+/*Headers: {{{1*/
+#ifdef HAVE_CONFIG_H
+	#include "config.h"
+#else
+#error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!"
+#endif
+
+#include <vector>
+#include <functional>
+#include <algorithm>
+#include <iostream>
+
+#include "./DataSet.h"
+#include "../shared/shared.h"
+#include "../include/include.h"
+#include "../EnumDefinitions/EnumDefinitions.h"
+
+using namespace std;
+/*}}}*/
+
+/*Object constructors and destructor*/
+/*FUNCTION Parameters::Parameters(){{{1*/
+Parameters::Parameters(){
+	return;
+}
+/*}}}*/
+/*FUNCTION Parameters::Parameters(int in_enum){{{1*/
+Parameters::Parameters(int in_enum): DataSet(in_enum){
+	//do nothing;
+	return;
+}
+/*}}}*/
+/*FUNCTION Parameters::~Parameters(){{{1*/
+Parameters::~Parameters(){
+	return;
+}
+/*}}}*/
+
+/*Object management*/
+/*FUNCTION Parameters::FindParam(bool* pbool,int enum_type){{{1*/
+int   Parameters::FindParam(bool* pbool,int enum_type){
+	
+	/*Go through a dataset, and find a Param* object 
+	 *which parameter name is "name" : */
+	
+	vector<Object*>::iterator object;
+	Param* param=NULL;
+
+	int found=0;
+
+	for ( object=objects.begin() ; object < objects.end(); object++ ){
+
+		/*Ok, this object is a parameter, recover it and ask which name it has: */
+		param=(Param*)(*object);
+
+		if(param->EnumType()==enum_type){
+			/*Ok, this is the one! Recover the value of this parameter: */
+			param->GetParameterValue(pbool);
+			found=1;
+			break;
+		}
+	}
+	return found;
+}
+/*}}}*/
+/*FUNCTION Parameters::FindParam(int* pinteger,int enum_type){{{1*/
+int   Parameters::FindParam(int* pinteger,int enum_type){
+	
+	/*Go through a dataset, and find a Param* object 
+	 *which parameter name is "name" : */
+	
+	vector<Object*>::iterator object;
+	Param* param=NULL;
+
+	int found=0;
+
+	for ( object=objects.begin() ; object < objects.end(); object++ ){
+
+		/*Ok, this object is a parameter, recover it and ask which name it has: */
+		param=(Param*)(*object);
+
+		if(param->EnumType()==enum_type){
+			/*Ok, this is the one! Recover the value of this parameter: */
+			param->GetParameterValue(pinteger);
+			found=1;
+			break;
+		}
+	}
+	return found;
+}
+/*}}}*/
+/*FUNCTION Parameters::FindParam(double* pscalar, int enum_type){{{1*/
+int   Parameters::FindParam(double* pscalar, int enum_type){
+	
+	/*Go through a dataset, and find a Param* object 
+	 *which parameter name is "name" : */
+	
+	vector<Object*>::iterator object;
+	Param* param=NULL;
+
+	int found=0;
+
+	for ( object=objects.begin() ; object < objects.end(); object++ ){
+
+		/*Ok, this object is a parameter, recover it and ask which name it has: */
+		param=(Param*)(*object);
+
+		if(param->EnumType()==enum_type){
+			/*Ok, this is the one! Recover the value of this parameter: */
+			param->GetParameterValue(pscalar);
+			found=1;
+			break;
+		}
+	}
+	return found;
+}
+/*}}}*/
+/*FUNCTION Parameters::FindParam(char** pstring,int enum_type){{{1*/
+int   Parameters::FindParam(char** pstring,int enum_type){
+	
+	/*Go through a dataset, and find a Param* object 
+	 *which parameter name is "name" : */
+	
+	vector<Object*>::iterator object;
+	Param* param=NULL;
+
+	int found=0;
+
+	for ( object=objects.begin() ; object < objects.end(); object++ ){
+
+		/*Ok, this object is a parameter, recover it and ask which name it has: */
+		param=(Param*)(*object);
+
+		if(param->EnumType()==enum_type){
+			/*Ok, this is the one! Recover the value of this parameter: */
+			param->GetParameterValue(pstring);
+			found=1;
+			break;
+		}
+	}
+	return found;
+
+}
+/*}}}*/
+/*FUNCTION Parameters::FindParam(char*** pstringarray,int* pM,int enum_type){{{1*/
+int   Parameters::FindParam(char*** pstringarray,int* pM,int enum_type){
+	
+	/*Go through a dataset, and find a Param* object 
+	 *which parameter name is "name" : */
+	
+	vector<Object*>::iterator object;
+	Param* param=NULL;
+
+	int found=0;
+
+	for ( object=objects.begin() ; object < objects.end(); object++ ){
+
+		/*Ok, this object is a parameter, recover it and ask which name it has: */
+		param=(Param*)(*object);
+
+		if(param->EnumType()==enum_type){
+			/*Ok, this is the one! Recover the value of this parameter: */
+			param->GetParameterValue(pstringarray,pM);
+			found=1;
+			break;
+		}
+	}
+	return found;
+
+}
+/*}}}*/
+/*FUNCTION Parameters::FindParam(double** pdoublearray,int* pM,int enum_type){{{1*/
+int   Parameters::FindParam(double** pdoublearray,int* pM, int enum_type){
+	
+	/*Go through a dataset, and find a Param* object 
+	 *which parameter name is "name" : */
+	
+	vector<Object*>::iterator object;
+	Param* param=NULL;
+
+	int found=0;
+
+	for ( object=objects.begin() ; object < objects.end(); object++ ){
+
+		/*Ok, this object is a parameter, recover it and ask which name it has: */
+		param=(Param*)(*object);
+
+		if(param->EnumType()==enum_type){
+			/*Ok, this is the one! Recover the value of this parameter: */
+			param->GetParameterValue(pdoublearray,pM);
+			found=1;
+			break;
+		}
+	}
+	return found;
+
+}
+/*}}}*/
+/*FUNCTION Parameters::FindParam(double** pdoublearray,int* pM, int* pN,int enum_type){{{1*/
+int   Parameters::FindParam(double** pdoublearray,int* pM, int* pN,int enum_type){
+	
+	/*Go through a dataset, and find a Param* object 
+	 *which parameter name is "name" : */
+	
+	vector<Object*>::iterator object;
+	Param* param=NULL;
+
+	int found=0;
+
+	for ( object=objects.begin() ; object < objects.end(); object++ ){
+
+		/*Ok, this object is a parameter, recover it and ask which name it has: */
+		param=(Param*)(*object);
+
+		if(param->EnumType()==enum_type){
+			/*Ok, this is the one! Recover the value of this parameter: */
+			param->GetParameterValue(pdoublearray,pM,pN);
+			found=1;
+			break;
+		}
+	}
+	return found;
+
+}
+/*}}}*/
+/*FUNCTION Parameters::FindParam(Vec* pvec,int enum_type){{{1*/
+int   Parameters::FindParam(Vec* pvec,int enum_type){
+	
+	/*Go through a dataset, and find a Param* object 
+	 *which parameter name is "name" : */
+	
+	vector<Object*>::iterator object;
+	Param* param=NULL;
+
+	int found=0;
+
+	for ( object=objects.begin() ; object < objects.end(); object++ ){
+
+		/*Ok, this object is a parameter, recover it and ask which name it has: */
+		param=(Param*)(*object);
+
+		if(param->EnumType()==enum_type){
+			/*Ok, this is the one! Recover the value of this parameter: */
+			param->GetParameterValue(pvec);
+			found=1;
+			break;
+		}
+	}
+	return found;
+
+}
+/*}}}*/
+/*FUNCTION Parameters::FindParam(Mat* pmat,int enum_type){{{1*/
+int   Parameters::FindParam(Mat* pmat,int enum_type){
+	
+	/*Go through a dataset, and find a Param* object 
+	 *which parameter name is "name" : */
+	
+	vector<Object*>::iterator object;
+	Param* param=NULL;
+
+	int found=0;
+
+	for ( object=objects.begin() ; object < objects.end(); object++ ){
+
+		/*Ok, this object is a parameter, recover it and ask which name it has: */
+		param=(Param*)(*object);
+
+		if(param->EnumType()==enum_type){
+			/*Ok, this is the one! Recover the value of this parameter: */
+			param->GetParameterValue(pmat);
+			found=1;
+			break;
+		}
+	}
+	return found;
+
+}
+/*}}}*/
+
+/*FUNCTION Parameters::SetParam(bool boolean,int enum_type);{{{1*/
+void   Parameters::SetParam(bool boolean,int enum_type){
+
+	Param* param=NULL;
+	
+	/*first, figure out if the param has already been created: */
+	param=(Param*)this->FindParamObject(enum_type);
+
+	if(param) param->SetValue(boolean); //already exists, just set it.
+	else this->AddObject(new BoolParam(enum_type,boolean)); //just add the new parameter.
+}
+/*}}}*/
+/*FUNCTION Parameters::SetParam(int integer,int enum_type);{{{1*/
+void   Parameters::SetParam(int integer,int enum_type){
+
+	Param* param=NULL;
+	
+	/*first, figure out if the param has already been created: */
+	param=(Param*)this->FindParamObject(enum_type);
+
+	if(param) param->SetValue(integer); //already exists, just set it.
+	else this->AddObject(new IntParam(enum_type,integer)); //just add the new parameter.
+}
+/*}}}*/
+/*FUNCTION Parameters::SetParam(double scalar,int enum_type);{{{1*/
+void   Parameters::SetParam(double scalar,int enum_type){
+
+	Param* param=NULL;
+	
+	/*first, figure out if the param has already been created: */
+	param=(Param*)this->FindParamObject(enum_type);
+
+	if(param) param->SetValue(scalar); //already exists, just set it.
+	else this->AddObject(new DoubleParam(enum_type,scalar)); //just add the new parameter.
+}
+/*}}}*/
+/*FUNCTION Parameters::SetParam(char* string,int enum_type);{{{1*/
+void   Parameters::SetParam(char* string,int enum_type){
+
+	Param* param=NULL;
+	
+	/*first, figure out if the param has already been created: */
+	param=(Param*)this->FindParamObject(enum_type);
+
+	if(param) param->SetValue(string); //already exists, just set it.
+	else this->AddObject(new StringParam(enum_type,string)); //just add the new parameter.
+}
+/*}}}*/
+/*FUNCTION Parameters::SetParam(char** stringarray,int M, int enum_type);{{{1*/
+void   Parameters::SetParam(char** stringarray,int M, int enum_type){
+
+	Param* param=NULL;
+	
+	/*first, figure out if the param has already been created: */
+	param=(Param*)this->FindParamObject(enum_type);
+
+	if(param) param->SetValue(stringarray,M); //already exists, just set it.
+	else this->AddObject(new StringArrayParam(enum_type,stringarray,M)); //just add the new parameter.
+}
+/*}}}*/
+/*FUNCTION Parameters::SetParam(double* doublearray,int M,int enum_type);{{{1*/
+void   Parameters::SetParam(double* doublearray,int M, int enum_type){
+
+	Param* param=NULL;
+	
+	/*first, figure out if the param has already been created: */
+	param=(Param*)this->FindParamObject(enum_type);
+
+	if(param) param->SetValue(doublearray,M); //already exists, just set it.
+	else this->AddObject(new DoubleVecParam(enum_type,doublearray,M)); //just add the new parameter.
+}
+/*}}}*/
+/*FUNCTION Parameters::SetParam(double* doublearray,int M,int N, int enum_type);{{{1*/
+void   Parameters::SetParam(double* doublearray,int M, int N, int enum_type){
+
+	Param* param=NULL;
+	
+	/*first, figure out if the param has already been created: */
+	param=(Param*)this->FindParamObject(enum_type);
+
+	if(param) param->SetValue(doublearray,M,N); //already exists, just set it.
+	else this->AddObject(new DoubleMatParam(enum_type,doublearray,M,N)); //just add the new parameter.
+}
+/*}}}*/
+/*FUNCTION Parameters::SetParam(Vec vector,int enum_type);{{{1*/
+void   Parameters::SetParam(Vec vector,int enum_type){
+
+	Param* param=NULL;
+	
+	/*first, figure out if the param has already been created: */
+	param=(Param*)this->FindParamObject(enum_type);
+
+	if(param) param->SetValue(vector); //already exists, just set it.
+	else this->AddObject(new PetscVecParam(enum_type,vector)); //just add the new parameter.
+}
+/*}}}*/
+/*FUNCTION Parameters::SetParam(Mat matrix,int enum_type);{{{1*/
+void   Parameters::SetParam(Mat matrix,int enum_type){
+
+	Param* param=NULL;
+	
+	/*first, figure out if the param has already been created: */
+	param=(Param*)this->FindParamObject(enum_type);
+
+	if(param) param->SetValue(matrix); //already exists, just set it.
+	else this->AddObject(new PetscMatParam(enum_type,matrix)); //just add the new parameter.
+}
+/*}}}*/
+
+/*FUNCTION Parameters::FindParamObject{{{1*/
+Object* Parameters::FindParamObject(int enum_type){
+
+	/*Go through a dataset, and find a Param* object 
+	 *which parameter name is "name" : */
+
+	vector<Object*>::iterator object;
+	Param* param=NULL;
+
+	for ( object=objects.begin() ; object < objects.end(); object++ ){
+
+		/*Ok, this object is a parameter, recover it and ask which name it has: */
+		param=(Param*)(*object);
+
+		if(param->EnumType()==enum_type){
+			/*Ok, this is the one! Return the object: */
+			return (*object);
+		}
+	}
+	return NULL;
+}
+/*}}}*/
Index: /issm/trunk/src/c/Container/Results.cpp
===================================================================
--- /issm/trunk/src/c/Container/Results.cpp	(revision 4235)
+++ /issm/trunk/src/c/Container/Results.cpp	(revision 4235)
@@ -0,0 +1,122 @@
+/*
+ * \file Results.c
+ * \brief: implementation of the Results class, derived from DataSet class
+ */
+
+/*Headers: {{{1*/
+#ifdef HAVE_CONFIG_H
+	#include "config.h"
+#else
+#error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!"
+#endif
+
+#include <vector>
+#include <functional>
+#include <algorithm>
+#include <iostream>
+
+#include "./DataSet.h"
+#include "../shared/shared.h"
+#include "../include/include.h"
+#include "../EnumDefinitions/EnumDefinitions.h"
+
+using namespace std;
+/*}}}*/
+
+/*Object constructors and destructor*/
+/*FUNCTION Results::Results(){{{1*/
+Results::Results(){
+	return;
+}
+/*}}}*/
+/*FUNCTION Results::Results(int in_enum){{{1*/
+Results::Results(int in_enum): DataSet(in_enum) {
+	//do nothing;
+	return;
+}
+/*}}}*/
+/*FUNCTION Results::~Results(){{{1*/
+Results::~Results(){
+	return;
+}
+/*}}}*/
+
+/*Object management*/
+/*FUNCTION Results::SpawnBeamResults{{{1*/
+Results* Results::SpawnBeamResults(int* indices){
+
+	/*Intermediary*/
+	vector<Object*>::iterator object;
+	ElementResult* resultin=NULL;
+	ElementResult* resultout=NULL;
+
+	/*Output*/
+	Results* newresults=new Results();
+
+	/*Go through results and call Spawn function*/
+	for ( object=objects.begin() ; object < objects.end(); object++ ){
+
+		/*Create new result*/
+		resultin=(ElementResult*)(*object); 
+		resultout=resultin->SpawnBeamElementResult(indices);
+
+		/*Add result to new results*/
+		newresults->AddObject((Object*)resultout);
+	}
+
+	/*Assign output pointer*/
+	return newresults;
+}
+/*}}}*/
+/*FUNCTION Results::SpawnSingResults{{{1*/
+Results* Results::SpawnSingResults(int index){
+
+	/*Intermediary*/
+	vector<Object*>::iterator object;
+	ElementResult* resultin=NULL;
+	ElementResult* resultout=NULL;
+
+	/*Output*/
+	Results* newresults=new Results();
+
+	/*Go through results and call Spawn function*/
+	for ( object=objects.begin() ; object < objects.end(); object++ ){
+
+		/*Create new result*/
+		resultin=(ElementResult*)(*object); 
+		resultout=resultin->SpawnSingElementResult(index);
+
+		/*Add result to new results*/
+		newresults->AddObject((Object*)resultout);
+	}
+
+	/*Assign output pointer*/
+	return newresults;
+}
+/*}}}*/
+/*FUNCTION Results::SpawnTriaResults{{{1*/
+Results* Results::SpawnTriaResults(int* indices){
+
+	/*Intermediary*/
+	vector<Object*>::iterator object;
+	ElementResult* resultin=NULL;
+	ElementResult* resultout=NULL;
+
+	/*Output*/
+	Results* newresults=new Results();
+
+	/*Go through results and call Spawn function*/
+	for ( object=objects.begin() ; object < objects.end(); object++ ){
+
+		/*Create new result*/
+		resultin=(ElementResult*)(*object); 
+		resultout=resultin->SpawnTriaElementResult(indices);
+
+		/*Add result to new results*/
+		newresults->AddObject((Object*)resultout);
+	}
+
+	/*Assign output pointer*/
+	return newresults;
+}
+/*}}}*/
Index: /issm/trunk/src/c/Container/Vertices.cpp
===================================================================
--- /issm/trunk/src/c/Container/Vertices.cpp	(revision 4235)
+++ /issm/trunk/src/c/Container/Vertices.cpp	(revision 4235)
@@ -0,0 +1,229 @@
+/*
+ * \file Vertices.c
+ * \brief: implementation of the Vertices class, derived from DataSet class
+ */
+
+/*Headers: {{{1*/
+#ifdef HAVE_CONFIG_H
+	#include "config.h"
+#else
+#error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!"
+#endif
+
+#include <vector>
+#include <functional>
+#include <algorithm>
+#include <iostream>
+
+#include "./DataSet.h"
+#include "../shared/shared.h"
+#include "../include/include.h"
+#include "../EnumDefinitions/EnumDefinitions.h"
+
+using namespace std;
+/*}}}*/
+
+/*Object constructors and destructor*/
+/*FUNCTION Vertices::Vertices(){{{1*/
+Vertices::Vertices(){
+	return;
+}
+/*}}}*/
+/*FUNCTION Vertices::Vertices(int in_enum){{{1*/
+Vertices::Vertices(int in_enum): DataSet(in_enum){
+	//do nothing;
+	return;
+}
+/*}}}*/
+/*FUNCTION Vertices::~Vertices(){{{1*/
+Vertices::~Vertices(){
+	return;
+}
+/*}}}*/
+
+/*Numerics management*/
+/*FUNCTION Vertices::CreatePartitioningVector{{{1*/
+void  Vertices::CreatePartitioningVector(Vec* ppartition,int numberofobjects){
+
+	int i;
+
+	/*output: */
+	Vec partition=NULL;
+
+	/*Create partition vector: */
+	partition=NewVec(numberofobjects);
+
+	/*Go through all objects, and ask each object to plug its doflist in 
+	 * partition: */
+
+	for(i=0;i<this->Size();i++){
+		Vertex* vertex=(Vertex*)this->GetObjectByOffset(i);
+		vertex->CreatePartition(partition);
+	}
+
+	/*Assemble the petsc vector: */
+	VecAssemblyBegin(partition);
+	VecAssemblyEnd(partition);
+
+	/*Assign output pointers: */
+	*ppartition=partition;
+}
+/*}}}*/
+/*FUNCTION Vertices::DistributeDofs{{{1*/
+void  Vertices::DistributeDofs(int numberofobjects,int numberofdofsperobject){
+
+	extern int num_procs;
+	extern int my_rank;
+
+	int  i;
+	
+	int  dofcount=0;
+	int* alldofcount=NULL;
+	int* truedofs=NULL;
+	int* alltruedofs=NULL;
+
+	/*Go through objects, and distribute dofs locally, from 0 to numberofdofsperobject: */
+	for (i=0;i<this->Size();i++){
+		Vertex* vertex=(Vertex*)this->GetObjectByOffset(i);
+		vertex->DistributeDofs(&dofcount);
+	}
+
+	/*Ok, now every object has distributed dofs, but locally, and with a dof count starting from 
+	 *0. This means the dofs between all the cpus are not synchronized! We need to synchronize all 
+	 *dof on all cpus, and use those to update the dofs of every object: */
+
+	alldofcount=(int*)xmalloc(num_procs*sizeof(int));
+	MPI_Gather(&dofcount,1,MPI_INT,alldofcount,1,MPI_INT,0,MPI_COMM_WORLD);
+	MPI_Bcast(alldofcount,num_procs,MPI_INT,0,MPI_COMM_WORLD);
+
+	/*Ok, now every cpu should start its own dof count at the end of the dofcount 
+	 * from cpu-1. : */
+	dofcount=0;
+	if(my_rank==0){
+		dofcount=0;
+	}
+	else{
+		for(i=0;i<my_rank;i++){
+			dofcount+=alldofcount[i];
+		}
+	}
+
+
+	/*Ok, now every cpu knows where his dofs should start. Update the dof count: */
+	for (i=0;i<this->Size();i++){
+		Vertex* vertex=(Vertex*)this->GetObjectByOffset(i);
+		vertex->OffsetDofs(dofcount);
+	}
+
+	/*Finally, remember that cpus may have skipped some objects, because they were clones. For every 
+	 * object that is not a clone, tell them to show their dofs, so that later on, they can get picked 
+	 * up by their clones: */
+	truedofs=(int*)xcalloc(numberofobjects*numberofdofsperobject,sizeof(int));
+	alltruedofs=(int*)xcalloc(numberofobjects*numberofdofsperobject,sizeof(int));
+
+	for (i=0;i<this->Size();i++){
+		Vertex* vertex=(Vertex*)this->GetObjectByOffset(i);
+		vertex->ShowTrueDofs(truedofs);
+	}
+	
+	MPI_Allreduce ( (void*)truedofs,(void*)alltruedofs,numberofobjects*numberofdofsperobject,MPI_INT,MPI_MAX,MPI_COMM_WORLD);
+
+	/*Ok, now every cpu knows the true dofs of everyone else that is not a clone. Let the clones recover those true dofs: */
+	for (i=0;i<this->Size();i++){
+		Vertex* vertex=(Vertex*)this->GetObjectByOffset(i);
+		vertex->UpdateCloneDofs(alltruedofs);
+	}
+
+	/* Free ressources: */
+	xfree((void**)&alldofcount);
+	xfree((void**)&truedofs);
+	xfree((void**)&alltruedofs);
+
+}
+/*}}}*/
+/*FUNCTION Vertices::FlagClones{{{1*/
+void  Vertices::FlagClones(int numberofobjects){
+
+	int i;
+	extern int num_procs;
+
+	int* ranks=NULL;
+	int* minranks=NULL;
+
+	/*Allocate ranks: */
+	ranks=(int*)xmalloc(numberofobjects*sizeof(int));
+	minranks=(int*)xmalloc(numberofobjects*sizeof(int));
+
+	for(i=0;i<numberofobjects;i++)ranks[i]=num_procs; //no cpu can have rank num_procs. This is the maximum limit.
+
+	/*Now go through all our objects and ask them to report to who they belong (which rank): */
+	Ranks(ranks);
+
+	/*We need to take the minimum rank for each vertex, and every cpu needs to get that result. That way, 
+	 * when we start building the dof list for all vertexs, a cpu can check whether its vertex already has been 
+	 * dealt with by another cpu. We take the minimum because we are going to manage dof assignment in increasing 
+	 * order of cpu rank. This is also why we initialized this array to num_procs.*/
+	MPI_Allreduce ( (void*)ranks,(void*)minranks,numberofobjects,MPI_INT,MPI_MIN,MPI_COMM_WORLD);
+
+	/*Now go through all objects, and use minranks to flag which objects are cloned: */
+	for(i=0;i<this->Size();i++){
+		/*For this object, decide whether it is a clone: */
+		Vertex* vertex=(Vertex*)this->GetObjectByOffset(i);
+		vertex->SetClone(minranks);
+	}
+
+	/*Free ressources: */
+	xfree((void**)&ranks); 
+	xfree((void**)&minranks);
+
+}
+/*}}}*/
+/*FUNCTION Vertices::NumberOfVertices{{{1*/
+int Vertices::NumberOfVertices(void){
+
+	int i;
+
+	int max_sid=0;
+	int sid;
+	int vertex_max_sid;
+
+	for(i=0;i<this->Size();i++){
+		
+		Vertex* vertex=(Vertex*)this->GetObjectByOffset(i);
+		sid=vertex->Sid();
+		if (sid>max_sid)max_sid=sid;
+	}
+
+	#ifdef _PARALLEL_
+	MPI_Reduce (&max_sid,&vertex_max_sid,1,MPI_INT,MPI_MAX,0,MPI_COMM_WORLD );
+	MPI_Bcast(&vertex_max_sid,1,MPI_INT,0,MPI_COMM_WORLD);
+	max_sid=vertex_max_sid;
+	#endif
+
+	/*sid starts at 0*/
+	max_sid++;
+
+	/*return:*/
+	return max_sid;
+}
+/*}}}*/
+/*FUNCTION Vertices::Ranks{{{1*/
+void   Vertices::Ranks(int* ranks){
+
+	/*Go through a dataset, and for each object, report it cpu: */
+
+	int i;
+	int rank;
+	int sid;
+	
+	for(i=0;i<this->Size();i++){
+		Vertex* vertex=(Vertex*)this->GetObjectByOffset(i);
+		rank=vertex->MyRank();
+		sid=vertex->Sid();
+		
+		/*Plug rank into ranks, according to id: */
+		ranks[sid]=rank; 
+	}
+	return;
+}
+/*}}}*/