Index: /issm/trunk/src/c/objects/Elements/Beam.cpp
===================================================================
--- /issm/trunk/src/c/objects/Elements/Beam.cpp	(revision 4010)
+++ /issm/trunk/src/c/objects/Elements/Beam.cpp	(revision 4011)
@@ -23,134 +23,4 @@
 	this->parameters=NULL;
 	return;
-}
-/*}}}*/
-/*FUNCTION Beam::Beam(int id, Hook* hnodes, Hook* hmatice, Hook* hmatpar, Parameters* beam_parameters, ElementProperties* properties){{{1*/
-Beam::Beam(int beam_id,Hook* beam_hnodes, Hook* beam_hmatice, Hook* beam_hmatpar, Parameters* beam_parameters, Inputs* beam_inputs):
-	hnodes(beam_hnodes),
-	hmatice(beam_hmatice),
-	hmatpar(beam_hmatpar)
-{
-
-	/*all the initialization has been done by the initializer, just fill in the id: */
-	this->id=beam_id;
-	if(beam_inputs){
-		this->inputs=(Inputs*)beam_inputs->Copy();
-	}
-	else{
-		this->inputs=new Inputs();
-	}
-	/*point parameters: */
-	this->parameters=beam_parameters;
-}
-/*}}}*/
-/*FUNCTION Beam::Beam(int id, int i, IoModel* iomodel){{{1*/
-Beam::Beam(int beam_id, int index,IoModel* iomodel){
-
-	int i;
-
-	/*beam constructor input: */
-	int   beam_matice_id;
-	int   beam_matpar_id;
-	int   beam_node_ids[2];
-	double nodeinputs[2];
-
-	/*id: */
-	this->id=beam_id;
-
-	/*hooks: */
-	ISSMASSERT(iomodel->uppernodes);
-	beam_matice_id=index+1; //refers to the corresponding material property card
-	beam_matpar_id=iomodel->numberofvertices2d*(iomodel->numlayers-1)+1;//refers to the corresponding matpar property card
-	beam_node_ids[0]=index+1;
-	beam_node_ids[1]=(int)iomodel->uppernodes[index]; //grid that lays right on top
-	
-	this->hnodes.Init(beam_node_ids,2);
-	this->hmatice.Init(&beam_matice_id,1);
-	this->hmatpar.Init(&beam_matpar_id,1);
-
-	//intialize inputs, and add as many inputs per element as requested: 
-	this->inputs=new Inputs();
-
-	if (iomodel->thickness) {
-		nodeinputs[0]=iomodel->thickness[index];
-		nodeinputs[1]=iomodel->thickness[(int)(iomodel->uppernodes[index]-1)];
-		this->inputs->AddInput(new BeamVertexInput(ThicknessEnum,nodeinputs));
-	}
-	if (iomodel->surface) {
-		nodeinputs[0]=iomodel->surface[index];
-		nodeinputs[1]=iomodel->surface[(int)(iomodel->uppernodes[index]-1)];
-		this->inputs->AddInput(new BeamVertexInput(SurfaceEnum,nodeinputs));
-	}	
-	if (iomodel->bed) {
-		nodeinputs[0]=iomodel->bed[index];
-		nodeinputs[1]=iomodel->bed[(int)(iomodel->uppernodes[index]-1)];
-		this->inputs->AddInput(new BeamVertexInput(BedEnum,nodeinputs));
-	}	
-	if (iomodel->drag_coefficient) {
-		nodeinputs[0]=iomodel->drag_coefficient[index];
-		nodeinputs[1]=iomodel->drag_coefficient[(int)(iomodel->uppernodes[index]-1)];
-		this->inputs->AddInput(new BeamVertexInput(DragCoefficientEnum,nodeinputs));
-	}
-	
-	if (iomodel->vx) {
-		nodeinputs[0]=iomodel->vx[index];
-		nodeinputs[1]=iomodel->vx[(int)(iomodel->uppernodes[index]-1)];
-		this->inputs->AddInput(new BeamVertexInput(VxEnum,nodeinputs));
-	}
-	if (iomodel->vy) {
-		nodeinputs[0]=iomodel->vy[index];
-		nodeinputs[1]=iomodel->vy[(int)(iomodel->uppernodes[index]-1)];
-		this->inputs->AddInput(new BeamVertexInput(VyEnum,nodeinputs));
-	}
-	if (iomodel->vz) {
-		nodeinputs[0]=iomodel->vz[index];
-		nodeinputs[1]=iomodel->vz[(int)(iomodel->uppernodes[index]-1)];
-		this->inputs->AddInput(new BeamVertexInput(VzEnum,nodeinputs));
-	}
-	if (iomodel->vx_obs) {
-		nodeinputs[0]=iomodel->vx_obs[index];
-		nodeinputs[1]=iomodel->vx_obs[(int)(iomodel->uppernodes[index]-1)];
-		this->inputs->AddInput(new BeamVertexInput(VxObsEnum,nodeinputs));
-	}
-	if (iomodel->vy_obs) {
-		nodeinputs[0]=iomodel->vy_obs[index];
-		nodeinputs[1]=iomodel->vy_obs[(int)(iomodel->uppernodes[index]-1)];
-		this->inputs->AddInput(new BeamVertexInput(VyObsEnum,nodeinputs));
-	}
-	if (iomodel->vz_obs) {
-		nodeinputs[0]=iomodel->vz_obs[index];
-		nodeinputs[1]=iomodel->vz_obs[(int)(iomodel->uppernodes[index]-1)];
-		this->inputs->AddInput(new BeamVertexInput(VzObsEnum,nodeinputs));
-	}
-	if (iomodel->weights) {
-		nodeinputs[0]=iomodel->weights[index];
-		nodeinputs[1]=iomodel->weights[(int)(iomodel->uppernodes[index]-1)];
-		this->inputs->AddInput(new BeamVertexInput(WeightsEnum,nodeinputs));
-	}	
-	if (iomodel->gridonbed) this->inputs->AddInput(new BoolInput(ElementOnBedEnum,(IssmBool)iomodel->gridonbed[index]));
-
-	/*default vx,vy and vz: */
-	if (!iomodel->vx && iomodel->vx_obs) {
-		nodeinputs[0]=iomodel->vx_obs[index];
-		nodeinputs[1]=iomodel->vx_obs[(int)(iomodel->uppernodes[index]-1)];
-		this->inputs->AddInput(new BeamVertexInput(VxEnum,nodeinputs));
-		this->inputs->AddInput(new BeamVertexInput(VxOldEnum,nodeinputs));
-	}
-	if (!iomodel->vy && iomodel->vy_obs) {
-		nodeinputs[0]=iomodel->vy_obs[index];
-		nodeinputs[1]=iomodel->vy_obs[(int)(iomodel->uppernodes[index]-1)];
-		this->inputs->AddInput(new BeamVertexInput(VyEnum,nodeinputs));
-		this->inputs->AddInput(new BeamVertexInput(VyOldEnum,nodeinputs));
-	}
-	if (!iomodel->vz && iomodel->vz_obs) {
-		nodeinputs[0]=iomodel->vz_obs[index];
-		nodeinputs[1]=iomodel->vz_obs[(int)(iomodel->uppernodes[index]-1)];
-		this->inputs->AddInput(new BeamVertexInput(VzEnum,nodeinputs));
-		this->inputs->AddInput(new BeamVertexInput(VzOldEnum,nodeinputs));
-	}
-
-	//this->parameters: we still can't point to it, it may not even exist. Configure will handle this.
-	this->parameters=NULL;
-
 }
 /*}}}*/
Index: /issm/trunk/src/c/objects/Elements/Beam.h
===================================================================
--- /issm/trunk/src/c/objects/Elements/Beam.h	(revision 4010)
+++ /issm/trunk/src/c/objects/Elements/Beam.h	(revision 4011)
@@ -36,6 +36,4 @@
 		/*constructors, destructors: {{{1*/
 		Beam();
-		Beam(int beam_id,Hook* beam_hnodes, Hook* beam_hmatice, Hook* beam_hmatpar, Parameters* beam_parameters, Inputs* beam_inputs);
-		Beam(int beam_id,int i, IoModel* iomodel);
 		~Beam();
 		/*}}}*/
Index: /issm/trunk/src/c/objects/Elements/Penta.cpp
===================================================================
--- /issm/trunk/src/c/objects/Elements/Penta.cpp	(revision 4010)
+++ /issm/trunk/src/c/objects/Elements/Penta.cpp	(revision 4011)
@@ -85,13 +85,26 @@
 Penta::Update(IoModel* iomodel,int analysis_counter){ 
 
+	/*Intermediaries*/
 	IssmInt i;
 	int penta_node_ids[6];
+	int penta_vertex_ids[6];
 	double nodeinputs[6];
 
+	/*Checks if debuging*/
+	/*{{{2*/
+	ISSMASSERT(iomodel->elements);
+	/*}}}*/
+
+	/*Recover vertices ids needed to initialize inputs*/
+	for(i=0;i<6;i++){ 
+		penta_vertex_ids[i]=(int)iomodel->elements[6*index+i]; //ids for vertices are in the elements array from Matlab
+	}
+
+	/*Recover nodes ids needed to initialize the node hook.*/
+	for(i=0;i<6;i++){ //go recover node ids, needed to initialize the node hook.
+		penta_node_ids[i]=(int)iomodel->elements[6*index+i]; //ids for vertices are in the elements array from Matlab
+	}
+
 	/*hooks: */
-	for(i=0;i<6;i++){ //go recover node ids, needed to initialize the node hook.
-		penta_node_ids[i]=(int)*(iomodel->elements+6*index+i); //ids for vertices are in the elements array from Matlab
-	}
-
 	this->SettHookNodes(penta_node_ids,analysis_counter); this->nodes=NULL; //set hook to nodes, for this analysis type
 
@@ -99,17 +112,17 @@
 	this->inputs=new Inputs();
 	if (iomodel->thickness) {
-		for(i=0;i<6;i++)nodeinputs[i]=iomodel->thickness[penta_node_ids[i]-1];
+		for(i=0;i<6;i++)nodeinputs[i]=iomodel->thickness[penta_vertex_ids[i]-1];
 		this->inputs->AddInput(new PentaVertexInput(ThicknessEnum,nodeinputs));
 	}
 	if (iomodel->surface) {
-		for(i=0;i<6;i++)nodeinputs[i]=iomodel->surface[penta_node_ids[i]-1];
+		for(i=0;i<6;i++)nodeinputs[i]=iomodel->surface[penta_vertex_ids[i]-1];
 		this->inputs->AddInput(new PentaVertexInput(SurfaceEnum,nodeinputs));
 	}
 	if (iomodel->bed) {
-		for(i=0;i<6;i++)nodeinputs[i]=iomodel->bed[penta_node_ids[i]-1];
+		for(i=0;i<6;i++)nodeinputs[i]=iomodel->bed[penta_vertex_ids[i]-1];
 		this->inputs->AddInput(new PentaVertexInput(BedEnum,nodeinputs));
 	}
 	if (iomodel->drag_coefficient) {
-		for(i=0;i<6;i++)nodeinputs[i]=iomodel->drag_coefficient[penta_node_ids[i]-1];
+		for(i=0;i<6;i++)nodeinputs[i]=iomodel->drag_coefficient[penta_vertex_ids[i]-1];
 		this->inputs->AddInput(new PentaVertexInput(DragCoefficientEnum,nodeinputs));
 
@@ -120,57 +133,57 @@
 	}
 	if (iomodel->melting_rate) {
-		for(i=0;i<6;i++)nodeinputs[i]=iomodel->melting_rate[penta_node_ids[i]-1]/iomodel->yts;
+		for(i=0;i<6;i++)nodeinputs[i]=iomodel->melting_rate[penta_vertex_ids[i]-1]/iomodel->yts;
 		this->inputs->AddInput(new PentaVertexInput(MeltingRateEnum,nodeinputs));
 	}
 	if (iomodel->accumulation_rate) {
-		for(i=0;i<6;i++)nodeinputs[i]=iomodel->accumulation_rate[penta_node_ids[i]-1]/iomodel->yts;
+		for(i=0;i<6;i++)nodeinputs[i]=iomodel->accumulation_rate[penta_vertex_ids[i]-1]/iomodel->yts;
 		this->inputs->AddInput(new PentaVertexInput(AccumulationRateEnum,nodeinputs));
 	}
 	if (iomodel->geothermalflux) {
-		for(i=0;i<6;i++)nodeinputs[i]=iomodel->geothermalflux[penta_node_ids[i]-1];
+		for(i=0;i<6;i++)nodeinputs[i]=iomodel->geothermalflux[penta_vertex_ids[i]-1];
 		this->inputs->AddInput(new PentaVertexInput(GeothermalFluxEnum,nodeinputs));
 	}	
 	if (iomodel->pressure) {
-		for(i=0;i<6;i++)nodeinputs[i]=iomodel->pressure[penta_node_ids[i]-1];
+		for(i=0;i<6;i++)nodeinputs[i]=iomodel->pressure[penta_vertex_ids[i]-1];
 		this->inputs->AddInput(new PentaVertexInput(PressureEnum,nodeinputs));
 	}
 	if (iomodel->temperature) {
-		for(i=0;i<6;i++)nodeinputs[i]=iomodel->temperature[penta_node_ids[i]-1];
+		for(i=0;i<6;i++)nodeinputs[i]=iomodel->temperature[penta_vertex_ids[i]-1];
 		this->inputs->AddInput(new PentaVertexInput(TemperatureEnum,nodeinputs));
 	}
 	if (iomodel->dhdt) {
-		for(i=0;i<6;i++)nodeinputs[i]=iomodel->dhdt[penta_node_ids[i]-1];
+		for(i=0;i<6;i++)nodeinputs[i]=iomodel->dhdt[penta_vertex_ids[i]-1];
 		this->inputs->AddInput(new PentaVertexInput(DhDtEnum,nodeinputs));
 	}
 	/*vx,vy and vz: */
 	if (iomodel->vx) {
-		for(i=0;i<6;i++)nodeinputs[i]=iomodel->vx[penta_node_ids[i]-1]/iomodel->yts;
+		for(i=0;i<6;i++)nodeinputs[i]=iomodel->vx[penta_vertex_ids[i]-1]/iomodel->yts;
 		this->inputs->AddInput(new PentaVertexInput(VxEnum,nodeinputs));
 		this->inputs->AddInput(new PentaVertexInput(VxOldEnum,nodeinputs));
 	}
 	if (iomodel->vy) {
-		for(i=0;i<6;i++)nodeinputs[i]=iomodel->vy[penta_node_ids[i]-1]/iomodel->yts;
+		for(i=0;i<6;i++)nodeinputs[i]=iomodel->vy[penta_vertex_ids[i]-1]/iomodel->yts;
 		this->inputs->AddInput(new PentaVertexInput(VyEnum,nodeinputs));
 		this->inputs->AddInput(new PentaVertexInput(VyOldEnum,nodeinputs));
 	}
 	if (iomodel->vz) {
-		for(i=0;i<6;i++)nodeinputs[i]=iomodel->vz[penta_node_ids[i]-1]/iomodel->yts;
+		for(i=0;i<6;i++)nodeinputs[i]=iomodel->vz[penta_vertex_ids[i]-1]/iomodel->yts;
 		this->inputs->AddInput(new PentaVertexInput(VzEnum,nodeinputs));
 		this->inputs->AddInput(new PentaVertexInput(VzOldEnum,nodeinputs));
 	}
 	if (iomodel->vx_obs) {
-		for(i=0;i<6;i++)nodeinputs[i]=iomodel->vx_obs[penta_node_ids[i]-1]/iomodel->yts;
+		for(i=0;i<6;i++)nodeinputs[i]=iomodel->vx_obs[penta_vertex_ids[i]-1]/iomodel->yts;
 		this->inputs->AddInput(new PentaVertexInput(VxObsEnum,nodeinputs));
 	}
 	if (iomodel->vy_obs) {
-		for(i=0;i<6;i++)nodeinputs[i]=iomodel->vy_obs[penta_node_ids[i]-1]/iomodel->yts;
+		for(i=0;i<6;i++)nodeinputs[i]=iomodel->vy_obs[penta_vertex_ids[i]-1]/iomodel->yts;
 		this->inputs->AddInput(new PentaVertexInput(VyObsEnum,nodeinputs));
 	}
 	if (iomodel->vz_obs) {
-		for(i=0;i<6;i++)nodeinputs[i]=iomodel->vz_obs[penta_node_ids[i]-1]/iomodel->yts;
+		for(i=0;i<6;i++)nodeinputs[i]=iomodel->vz_obs[penta_vertex_ids[i]-1]/iomodel->yts;
 		this->inputs->AddInput(new PentaVertexInput(VzObsEnum,nodeinputs));
 	}
 	if (iomodel->weights) {
-		for(i=0;i<6;i++)nodeinputs[i]=iomodel->weights[penta_node_ids[i]-1];
+		for(i=0;i<6;i++)nodeinputs[i]=iomodel->weights[penta_vertex_ids[i]-1];
 		this->inputs->AddInput(new PentaVertexInput(WeightsEnum,nodeinputs));
 	}
@@ -178,5 +191,5 @@
 	/*default vx,vy and vz: either observation or 0 */
 	if (!iomodel->vx){
-		if (iomodel->vx_obs) for(i=0;i<6;i++)nodeinputs[i]=iomodel->vx_obs[penta_node_ids[i]-1]/iomodel->yts;
+		if (iomodel->vx_obs) for(i=0;i<6;i++)nodeinputs[i]=iomodel->vx_obs[penta_vertex_ids[i]-1]/iomodel->yts;
 		else                 for(i=0;i<6;i++)nodeinputs[i]=0;
 		this->inputs->AddInput(new PentaVertexInput(VxEnum,nodeinputs));
@@ -184,5 +197,5 @@
 	}
 	if (!iomodel->vy){
-		if (iomodel->vy_obs) for(i=0;i<6;i++)nodeinputs[i]=iomodel->vy_obs[penta_node_ids[i]-1]/iomodel->yts;
+		if (iomodel->vy_obs) for(i=0;i<6;i++)nodeinputs[i]=iomodel->vy_obs[penta_vertex_ids[i]-1]/iomodel->yts;
 		else                 for(i=0;i<6;i++)nodeinputs[i]=0;
 		this->inputs->AddInput(new PentaVertexInput(VyEnum,nodeinputs));
@@ -190,5 +203,5 @@
 	}
 	if (!iomodel->vz){
-		if (iomodel->vz_obs) for(i=0;i<6;i++)nodeinputs[i]=iomodel->vz_obs[penta_node_ids[i]-1]/iomodel->yts;
+		if (iomodel->vz_obs) for(i=0;i<6;i++)nodeinputs[i]=iomodel->vz_obs[penta_vertex_ids[i]-1]/iomodel->yts;
 		else                 for(i=0;i<6;i++)nodeinputs[i]=0;
 		this->inputs->AddInput(new PentaVertexInput(VzEnum,nodeinputs));
Index: /issm/trunk/src/c/objects/Elements/Sing.cpp
===================================================================
--- /issm/trunk/src/c/objects/Elements/Sing.cpp	(revision 4010)
+++ /issm/trunk/src/c/objects/Elements/Sing.cpp	(revision 4011)
@@ -23,75 +23,4 @@
 	this->parameters=NULL;
 	return;
-}
-/*}}}*/
-/*FUNCTION Sing::Sing(int sing_id,Hook* sing_hnodes, Hook* sing_hmatice, Hook* sing_hmatpar, Hook* sing_hnumpar, Inputs* sing_inputs) {{{1*/
-Sing::Sing(int sing_id,Hook* sing_hnodes, Hook* sing_hmatice, Hook* sing_hmatpar, Parameters* sing_parameters,Inputs* sing_inputs):
-	hnodes(sing_hnodes),
-	hmatice(sing_hmatice),
-	hmatpar(sing_hmatpar)
-{
-
-	/*all the initialization has been done by the initializer, just fill in the id: */
-	this->id=sing_id;
-	if(sing_inputs){
-		this->inputs=(Inputs*)sing_inputs->Copy();
-	}
-	else{
-		this->inputs=new Inputs();
-	}
-	/*point parameters: */
-	this->parameters=sing_parameters;
-}
-/*}}}*/
-/*FUNCTION Sing::Sing(int sing_id, int i, IoModel* iomodel) {{{1*/
-Sing::Sing(int sing_id, int i, IoModel* iomodel){
-
-	int sing_matice_id;
-	int sing_matpar_id;
-	
-	int sing_g;
-	double sing_h;
-	double sing_k;
-
-	/*id: */
-	this->id=sing_id;
-
-	/*hooks: */
-	sing_matice_id=i+1;                        //refers to the corresponding material property card
-	sing_matpar_id=iomodel->numberofvertices+1;//refers to the corresponding matpar property card
-	sing_g=i+1;
-
-	this->hnodes.Init(&sing_g,1);
-	this->hmatice.Init(&sing_matice_id,1);
-	this->hmatpar.Init(&sing_matpar_id,1);
-
-	//intialize inputs, and add as many inputs per element as requested: 
-	this->inputs=new Inputs();
-
-	if (iomodel->thickness) this->inputs->AddInput(new SingVertexInput(ThicknessEnum,iomodel->thickness[i]));
-	if (iomodel->drag_coefficient) this->inputs->AddInput(new SingVertexInput(DragCoefficientEnum,iomodel->drag_coefficient[i]));
-	if (iomodel->vx) this->inputs->AddInput(new SingVertexInput(VxEnum,iomodel->vx[i]));
-	if (iomodel->vy) this->inputs->AddInput(new SingVertexInput(VyEnum,iomodel->vy[i]));
-	if (iomodel->vz) this->inputs->AddInput(new SingVertexInput(VzEnum,iomodel->vz[i]));
-	if (iomodel->vx_obs) this->inputs->AddInput(new SingVertexInput(VxObsEnum,iomodel->vx_obs[i]));
-	if (iomodel->vy_obs) this->inputs->AddInput(new SingVertexInput(VyObsEnum,iomodel->vy_obs[i]));
-	if (iomodel->vz_obs) this->inputs->AddInput(new SingVertexInput(VzObsEnum,iomodel->vz_obs[i]));
-	if (iomodel->weights) this->inputs->AddInput(new SingVertexInput(WeightsEnum,iomodel->weights[i]));
-	if (!iomodel->vx && iomodel->vx_obs){
-		this->inputs->AddInput(new SingVertexInput(VxEnum,iomodel->vx_obs[i]));
-		this->inputs->AddInput(new SingVertexInput(VxOldEnum,iomodel->vx_obs[i]));
-	}
-	if (!iomodel->vy && iomodel->vy_obs){
-		this->inputs->AddInput(new SingVertexInput(VyEnum,iomodel->vy_obs[i]));
-		this->inputs->AddInput(new SingVertexInput(VyOldEnum,iomodel->vy_obs[i]));
-	}
-	if (!iomodel->vz && iomodel->vz_obs){
-		this->inputs->AddInput(new SingVertexInput(VzEnum,iomodel->vz_obs[i]));
-		this->inputs->AddInput(new SingVertexInput(VzOldEnum,iomodel->vz_obs[i]));
-	}
-
-	//this->parameters: we still can't point to it, it may not even exist. Configure will handle this.
-	this->parameters=NULL;
-
 }
 /*}}}*/
Index: /issm/trunk/src/c/objects/Elements/Sing.h
===================================================================
--- /issm/trunk/src/c/objects/Elements/Sing.h	(revision 4010)
+++ /issm/trunk/src/c/objects/Elements/Sing.h	(revision 4011)
@@ -36,6 +36,4 @@
 		/*constructors, destructors: {{{1*/
 		Sing();
-		Sing(int sing_id,Hook* sing_hnodes, Hook* sing_hmatice, Hook* sing_hmatpar, Parameters* sing_parameters,Inputs* sing_inputs);
-		Sing(int sing_id, int i, IoModel* iomodel);
 		~Sing();
 		/*}}}*/
Index: /issm/trunk/src/c/objects/Elements/Tria.cpp
===================================================================
--- /issm/trunk/src/c/objects/Elements/Tria.cpp	(revision 4010)
+++ /issm/trunk/src/c/objects/Elements/Tria.cpp	(revision 4011)
@@ -68,43 +68,54 @@
 Tria::Update(IoModel* iomodel,int analysis_counter,int analysis_type){ //i is the element index
 
+	/*Intermediaries*/
 	int    i;
 	int    tria_node_ids[3];
+	int    tria_vertex_ids[3];
 	double nodeinputs[3];
 
-	/*hooks: */
-	//go recover node ids, needed to initialize the node hook.
+	/*Checks if debuging*/
+	/*{{{2*/
+	ISSMASSERT(iomodel->elements);
+	/*}}}*/
+
+	/*Recover vertices ids needed to initialize inputs*/
+	for(i=0;i<3;i++){ 
+		tria_vertex_ids[i]=(int)iomodel->elements[3*index+i]; //ids for vertices are in the elements array from Matlab
+	}
+
+	/*Recover nodes ids needed to initialize the node hook.*/
 	if (iomodel->analysis_type==Prognostic2AnalysisEnum || iomodel->analysis_type==Balancedthickness2AnalysisEnum){
 		/*Discontinuous Galerkin*/
-		tria_node_ids[0]=3*index+1;
-		tria_node_ids[1]=3*index+2;
-		tria_node_ids[2]=3*index+3;
+		tria_node_ids[0]=iomodel->nodecounter+3*index+1;
+		tria_node_ids[1]=iomodel->nodecounter+3*index+2;
+		tria_node_ids[2]=iomodel->nodecounter+3*index+3;
 	}
 	else{
 		/*Continuous Galerkin*/
-		ISSMASSERT(iomodel->elements);
 		for(i=0;i<3;i++){ 
-			tria_node_ids[i]=(int)*(iomodel->elements+3*index+i); //ids for vertices are in the elements array from Matlab
+			tria_node_ids[i]=iomodel->nodecounter+(int)*(iomodel->elements+3*index+i); //ids for vertices are in the elements array from Matlab
 		}
 	}
-	
+
+	/*hooks: */
 	this->SetHookNodes(tria_node_ids,analysis_counter); this->nodes=NULL; //set hook to nodes, for this analysis type
-
-	//intialize inputs, and add as many inputs per element as requested: 
+	
+	/*intialize inputs, and add as many inputs per element as requested:*/
 	this->inputs=new Inputs();
 	
 	if (iomodel->thickness) {
-		for(i=0;i<3;i++)nodeinputs[i]=iomodel->thickness[tria_node_ids[i]-1];
+		for(i=0;i<3;i++)nodeinputs[i]=iomodel->thickness[tria_vertex_ids[i]-1];
 		this->inputs->AddInput(new TriaVertexInput(ThicknessEnum,nodeinputs));
 	}
 	if (iomodel->surface) {
-		for(i=0;i<3;i++)nodeinputs[i]=iomodel->surface[tria_node_ids[i]-1];
+		for(i=0;i<3;i++)nodeinputs[i]=iomodel->surface[tria_vertex_ids[i]-1];
 		this->inputs->AddInput(new TriaVertexInput(SurfaceEnum,nodeinputs));
 	}
 	if (iomodel->bed) {
-		for(i=0;i<3;i++)nodeinputs[i]=iomodel->bed[tria_node_ids[i]-1];
+		for(i=0;i<3;i++)nodeinputs[i]=iomodel->bed[tria_vertex_ids[i]-1];
 		this->inputs->AddInput(new TriaVertexInput(BedEnum,nodeinputs));
 	}
 	if (iomodel->drag_coefficient) {
-		for(i=0;i<3;i++)nodeinputs[i]=iomodel->drag_coefficient[tria_node_ids[i]-1];
+		for(i=0;i<3;i++)nodeinputs[i]=iomodel->drag_coefficient[tria_vertex_ids[i]-1];
 		this->inputs->AddInput(new TriaVertexInput(DragCoefficientEnum,nodeinputs));
 
@@ -114,57 +125,57 @@
 	}
 	if (iomodel->melting_rate) {
-		for(i=0;i<3;i++)nodeinputs[i]=iomodel->melting_rate[tria_node_ids[i]-1]/iomodel->yts;
+		for(i=0;i<3;i++)nodeinputs[i]=iomodel->melting_rate[tria_vertex_ids[i]-1]/iomodel->yts;
 		this->inputs->AddInput(new TriaVertexInput(MeltingRateEnum,nodeinputs));
 	}
 	if (iomodel->accumulation_rate) {
-		for(i=0;i<3;i++)nodeinputs[i]=iomodel->accumulation_rate[tria_node_ids[i]-1]/iomodel->yts;
+		for(i=0;i<3;i++)nodeinputs[i]=iomodel->accumulation_rate[tria_vertex_ids[i]-1]/iomodel->yts;
 		this->inputs->AddInput(new TriaVertexInput(AccumulationRateEnum,nodeinputs));
 	}
 	if (iomodel->geothermalflux) {
-		for(i=0;i<3;i++)nodeinputs[i]=iomodel->geothermalflux[tria_node_ids[i]-1];
+		for(i=0;i<3;i++)nodeinputs[i]=iomodel->geothermalflux[tria_vertex_ids[i]-1];
 		this->inputs->AddInput(new TriaVertexInput(GeothermalFluxEnum,nodeinputs));
 	}
 	if (iomodel->dhdt) {
-		for(i=0;i<3;i++)nodeinputs[i]=iomodel->dhdt[tria_node_ids[i]-1];
+		for(i=0;i<3;i++)nodeinputs[i]=iomodel->dhdt[tria_vertex_ids[i]-1];
 		this->inputs->AddInput(new TriaVertexInput(DhDtEnum,nodeinputs));
 	}
 	if (iomodel->pressure) {
-		for(i=0;i<3;i++)nodeinputs[i]=iomodel->pressure[tria_node_ids[i]-1];
+		for(i=0;i<3;i++)nodeinputs[i]=iomodel->pressure[tria_vertex_ids[i]-1];
 		this->inputs->AddInput(new TriaVertexInput(PressureEnum,nodeinputs));
 	}
 	if (iomodel->temperature) {
-		for(i=0;i<3;i++)nodeinputs[i]=iomodel->temperature[tria_node_ids[i]-1];
+		for(i=0;i<3;i++)nodeinputs[i]=iomodel->temperature[tria_vertex_ids[i]-1];
 		this->inputs->AddInput(new TriaVertexInput(TemperatureEnum,nodeinputs));
 	}
 	/*vx,vy and vz: */
 	if (iomodel->vx) {
-		for(i=0;i<3;i++)nodeinputs[i]=iomodel->vx[tria_node_ids[i]-1]/iomodel->yts;
+		for(i=0;i<3;i++)nodeinputs[i]=iomodel->vx[tria_vertex_ids[i]-1]/iomodel->yts;
 		this->inputs->AddInput(new TriaVertexInput(VxEnum,nodeinputs));
 		this->inputs->AddInput(new TriaVertexInput(VxOldEnum,nodeinputs));
 	}
 	if (iomodel->vy) {
-		for(i=0;i<3;i++)nodeinputs[i]=iomodel->vy[tria_node_ids[i]-1]/iomodel->yts;
+		for(i=0;i<3;i++)nodeinputs[i]=iomodel->vy[tria_vertex_ids[i]-1]/iomodel->yts;
 		this->inputs->AddInput(new TriaVertexInput(VyEnum,nodeinputs));
 		this->inputs->AddInput(new TriaVertexInput(VyOldEnum,nodeinputs));
 	}
 	if (iomodel->vz) {
-		for(i=0;i<3;i++)nodeinputs[i]=iomodel->vz[tria_node_ids[i]-1]/iomodel->yts;
+		for(i=0;i<3;i++)nodeinputs[i]=iomodel->vz[tria_vertex_ids[i]-1]/iomodel->yts;
 		this->inputs->AddInput(new TriaVertexInput(VzEnum,nodeinputs));
 		this->inputs->AddInput(new TriaVertexInput(VzOldEnum,nodeinputs));
 	}
 	if (iomodel->vx_obs) {
-		for(i=0;i<3;i++)nodeinputs[i]=iomodel->vx_obs[tria_node_ids[i]-1]/iomodel->yts;
+		for(i=0;i<3;i++)nodeinputs[i]=iomodel->vx_obs[tria_vertex_ids[i]-1]/iomodel->yts;
 		this->inputs->AddInput(new TriaVertexInput(VxObsEnum,nodeinputs));
 	}
 	if (iomodel->vy_obs) {
-		for(i=0;i<3;i++)nodeinputs[i]=iomodel->vy_obs[tria_node_ids[i]-1]/iomodel->yts;
+		for(i=0;i<3;i++)nodeinputs[i]=iomodel->vy_obs[tria_vertex_ids[i]-1]/iomodel->yts;
 		this->inputs->AddInput(new TriaVertexInput(VyObsEnum,nodeinputs));
 	}
 	if (iomodel->vz_obs) {
-		for(i=0;i<3;i++)nodeinputs[i]=iomodel->vz_obs[tria_node_ids[i]-1]/iomodel->yts;
+		for(i=0;i<3;i++)nodeinputs[i]=iomodel->vz_obs[tria_vertex_ids[i]-1]/iomodel->yts;
 		this->inputs->AddInput(new TriaVertexInput(VzObsEnum,nodeinputs));
 	}
 	if (iomodel->weights) {
-		for(i=0;i<3;i++)nodeinputs[i]=iomodel->weights[tria_node_ids[i]-1];
+		for(i=0;i<3;i++)nodeinputs[i]=iomodel->weights[tria_vertex_ids[i]-1];
 		this->inputs->AddInput(new TriaVertexInput(WeightsEnum,nodeinputs));
 	}
@@ -172,5 +183,5 @@
 	/*default vx,vy and vz: either observation or 0 */
 	if(!iomodel->vx){
-		if (iomodel->vx_obs) for(i=0;i<3;i++)nodeinputs[i]=iomodel->vx_obs[tria_node_ids[i]-1]/iomodel->yts;
+		if (iomodel->vx_obs) for(i=0;i<3;i++)nodeinputs[i]=iomodel->vx_obs[tria_vertex_ids[i]-1]/iomodel->yts;
 		else                 for(i=0;i<3;i++)nodeinputs[i]=0;
 		this->inputs->AddInput(new TriaVertexInput(VxEnum,nodeinputs));
@@ -178,5 +189,5 @@
 	}
 	if(!iomodel->vy){
-		if (iomodel->vy_obs) for(i=0;i<3;i++)nodeinputs[i]=iomodel->vy_obs[tria_node_ids[i]-1]/iomodel->yts;
+		if (iomodel->vy_obs) for(i=0;i<3;i++)nodeinputs[i]=iomodel->vy_obs[tria_vertex_ids[i]-1]/iomodel->yts;
 		else                 for(i=0;i<3;i++)nodeinputs[i]=0;
 		this->inputs->AddInput(new TriaVertexInput(VyEnum,nodeinputs));
@@ -184,5 +195,5 @@
 	}
 	if(!iomodel->vz){
-		if (iomodel->vz_obs) for(i=0;i<3;i++)nodeinputs[i]=iomodel->vz_obs[tria_node_ids[i]-1]/iomodel->yts;
+		if (iomodel->vz_obs) for(i=0;i<3;i++)nodeinputs[i]=iomodel->vz_obs[tria_vertex_ids[i]-1]/iomodel->yts;
 		else                 for(i=0;i<3;i++)nodeinputs[i]=0;
 		this->inputs->AddInput(new TriaVertexInput(VzEnum,nodeinputs));