Index: /issm/trunk-jpl/src/c/Makefile.am
===================================================================
--- /issm/trunk-jpl/src/c/Makefile.am	(revision 13219)
+++ /issm/trunk-jpl/src/c/Makefile.am	(revision 13220)
@@ -46,5 +46,4 @@
 					./classes/objects/Load.h\
 					./classes/objects/Contour.h\
-					./classes/objects/Contour.cpp\
 					./classes/objects/OptArgs.h\
 					./classes/objects/OptPars.h\
@@ -62,5 +61,4 @@
 					./classes/objects/Node.cpp\
 					./classes/objects/Segment.h\
-					./classes/objects/Segment.cpp\
 					./classes/objects/Vertex.h\
 					./classes/objects/Vertex.cpp\
@@ -311,5 +309,4 @@
 					./modules/UpdateConstraintsx/UpdateConstraintsx.h\
 					./modules/UpdateConstraintsx/UpdateConstraintsx.cpp\
-					./modules/MeshPartitionx/MeshPartitionx.cpp\
 					./modules/MeshPartitionx/MeshPartitionx.h\
 					./modules/Reducevectorgtofx/Reducevectorgtofx.cpp\
@@ -844,5 +841,4 @@
 			./shared/Exp/IsInPoly.cpp\
 			./shared/Exp/IsInPolySerial.cpp\
-			./shared/Exp/DomainOutlineRead.cpp\
 			./shared/Exp/DomainOutlineWrite.cpp\
 			./shared/TriMesh/trimesh.h\
@@ -1014,6 +1010,10 @@
 libISSMModules_a_SOURCES = $(module_sources)
 libISSMModules_a_SOURCES += $(bamg_sources)
+if KRIGING
 libISSMModules_a_SOURCES += $(kriging_sources)
+endif
+if KML
 libISSMModules_a_SOURCES += $(kml_sources)
+endif
 libISSMModules_a_CXXFLAGS = $(ALLCXXFLAGS)
 libISSMModules_a_LIBADD = ./libISSMCore.a
Index: /issm/trunk-jpl/src/c/classes/matrix/Matrix.h
===================================================================
--- /issm/trunk-jpl/src/c/classes/matrix/Matrix.h	(revision 13219)
+++ /issm/trunk-jpl/src/c/classes/matrix/Matrix.h	(revision 13220)
@@ -78,7 +78,7 @@
 		/*FUNCTION Matrix(int M,int N,IssmDouble sparsity,int in_type){{{*/
 		#ifdef _HAVE_PETSC_
-		Matrix(int M,int N,IssmDouble sparsity,int in_type=PetscMatType){
+		Matrix(int M,int N,double sparsity,int in_type=PetscMatType){
 		#else
-		Matrix(int M,int N,IssmDouble sparsity,int in_type=SeqMatType){
+		Matrix(int M,int N,double sparsity,int in_type=SeqMatType){
 		#endif
 
@@ -105,7 +105,7 @@
 		/*FUNCTION Matrix(IssmDouble* serial_mat, int M,int N,IssmDouble sparsity,int in_type){{{*/
 		#ifdef _HAVE_PETSC_
-		Matrix(IssmDouble* serial_mat,int M,int N,IssmDouble sparsity,int in_type=PetscMatType){
+		Matrix(IssmPDouble* serial_mat,int M,int N,IssmPDouble sparsity,int in_type=PetscMatType){
 		#else
-		Matrix(IssmDouble* serial_mat,int M,int N,IssmDouble sparsity,int in_type=SeqMatType){
+		Matrix(IssmPDouble* serial_mat,int M,int N,IssmPDouble sparsity,int in_type=SeqMatType){
 		#endif
 
Index: sm/trunk-jpl/src/c/classes/objects/Contour.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/objects/Contour.cpp	(revision 13219)
+++ 	(revision )
@@ -1,95 +1,0 @@
-/*! \file Contour.c
- *  \sa Contour.h
- */
-
-#ifdef HAVE_CONFIG_H
-	#include <config.h>
-#else
-#error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!"
-#endif
-
-#include <string.h>
-#include "./objects.h"
-#include "../../include/include.h"
-#include "../../io/io.h"
-
-/*Contour constructors and destructors:*/
-/*FUNCTION Contour::Contour() default constructor {{{*/
-Contour::Contour(){
-	this->id=0;
-	this->nods=0;
-	this->x=NULL;
-	this->y=NULL;
-	this->closed=false;
-}
-/*}}}*/
-/*FUNCTION Contour::Contour(int pid, int nods, IssmDouble* x, IssmDouble* y,bool closed) {{{*/
-Contour::Contour(int pid,int pnods, IssmDouble* px, IssmDouble* py,bool pclosed){
-	
-	this->id=pid;
-	this->nods=pnods;
-	this->closed=pclosed;
-	if(nods){
-		this->x=xNew<IssmDouble>(nods);
-		xMemCpy<IssmDouble>(this->x,px,nods);
-		this->y=xNew<IssmDouble>(nods);
-		xMemCpy<IssmDouble>(this->y,py,nods);
-	}
-}
-/*}}}*/
-/*FUNCTION Contour::Contour() default constructor {{{*/
-Contour::~Contour(){
-	xDelete<IssmDouble>(this->x);
-	xDelete<IssmDouble>(this->y);
-}
-/*}}}*/
-
-
-/*Object virtual function resolutoin: */
-/*FUNCTION Contour::Echo(){{{*/
-void Contour::Echo(void){
-
-	int i;
-
-	_printLine_("Contour: " << id);
-	_printLine_("   nods: " << nods);
-	_printLine_("   closed: " << (closed?"true":"false"));
-	if(nods){
-	        _printLine_("   x,y:");
-		for(i=0;i<nods;i++){
-		        _printLine_(i << ": " << x[i] << "|" << y[i]);
-		}
-	}
-}
-/*}}}*/
-/*FUNCTION Contour::DeepEcho(){{{*/
-void Contour::DeepEcho(void){
-	this->Echo();
-}
-/*}}}*/
-/*FUNCTION Contour::Id(){{{*/
-int Contour::Id(void){
-	return id;
-}
-/*}}}*/
-/*FUNCTION Contour::MyRank{{{*/
-int    Contour::MyRank(void){ 
-	extern int my_rank;
-
-	return my_rank; 
-}
-/*}}}*/
-/*FUNCTION Contour::ObjectEnum{{{*/
-int Contour::ObjectEnum(void){
-
-	return ContourEnum;
-
-}
-/*}}}*/
-/*FUNCTION Contour::copy {{{*/
-Object* Contour::copy() {
-
-	return new Contour(*this); 
-
-}
-/*}}}*/
Index: /issm/trunk-jpl/src/c/classes/objects/Contour.h
===================================================================
--- /issm/trunk-jpl/src/c/classes/objects/Contour.h	(revision 13219)
+++ /issm/trunk-jpl/src/c/classes/objects/Contour.h	(revision 13220)
@@ -8,10 +8,14 @@
 /*Headers:*/
 /*{{{*/
+#include "../../include/include.h"
+#include "../../shared/Exceptions/exceptions.h"
+#include "../../shared/Alloc/xNewDelete.h"
+#include "../../shared/MemOps/xMemCpy.h"
+#include "../../io/io.h"
+#include "../../EnumDefinitions/EnumDefinitions.h"
 #include "./Object.h"
-#include "../../shared/Exceptions/exceptions.h"
-#include "../../toolkits/toolkits.h"
-#include "../../include/include.h"
 /*}}}*/
 
+template <class doubletype>
 class Contour: public Object{
 
@@ -20,23 +24,88 @@
 		int     id;
 		int	  nods;  //number of vertices in the contour
-		IssmDouble* x;
-		IssmDouble* y;
+		doubletype* x;
+		doubletype* y;
 		bool    closed; //is this contour closed?
 
-		/*Contour constructors, destructors {{{*/
-		Contour();
-		Contour(int id, int nods, IssmDouble* x, IssmDouble* y,bool closed);
-		~Contour();
+		/*Contour constructors, destructors :*/
+		/*FUNCTION Contour() default constructor {{{*/
+		Contour(){
+			this->id=0;
+			this->nods=0;
+			this->x=NULL;
+			this->y=NULL;
+			this->closed=false;
+		}
 		/*}}}*/
-		/*Object virtual functions{{{*/
-		void  Echo(void);
-		void  DeepEcho(void);
-		int   Id(void);
-		int   MyRank(void);
-		int   ObjectEnum(void);
-		Object* copy(void);
+		/*FUNCTION Contour(int pid, int nods, doubletype* x, doubletype* y,bool closed) {{{*/
+		Contour(int pid,int pnods, doubletype* px, doubletype* py,bool pclosed){
+
+			this->id=pid;
+			this->nods=pnods;
+			this->closed=pclosed;
+			if(nods){
+				this->x=xNew<doubletype>(nods);
+				xMemCpy<doubletype>(this->x,px,nods);
+				this->y=xNew<doubletype>(nods);
+				xMemCpy<doubletype>(this->y,py,nods);
+			}
+		}
+		/*}}}*/
+		/*FUNCTION Contour() default constructor {{{*/
+		~Contour(){
+			xDelete<doubletype>(this->x);
+			xDelete<doubletype>(this->y);
+		}
 		/*}}}*/
 
 
+		/*Object virtual function resolutoin: */
+		/*FUNCTION Echo(){{{*/
+		void Echo(void){
+
+			int i;
+
+			_printLine_("Contour: " << id);
+			_printLine_("   nods: " << nods);
+			_printLine_("   closed: " << (closed?"true":"false"));
+			if(nods){
+				_printLine_("   x,y:");
+				for(i=0;i<nods;i++){
+					_printLine_(i << ": " << x[i] << "|" << y[i]);
+				}
+			}
+		}
+		/*}}}*/
+		/*FUNCTION DeepEcho(){{{*/
+		void DeepEcho(void){
+			this->Echo();
+		}
+		/*}}}*/
+		/*FUNCTION Id(){{{*/
+		int Id(void){
+			return id;
+		}
+		/*}}}*/
+		/*FUNCTION MyRank{{{*/
+		int    MyRank(void){ 
+			extern int my_rank;
+
+			return my_rank; 
+		}
+		/*}}}*/
+		/*FUNCTION ObjectEnum{{{*/
+		int ObjectEnum(void){
+
+			return ContourEnum;
+
+		}
+		/*}}}*/
+		/*FUNCTION copy {{{*/
+		Object* copy() {
+
+			return new Contour(*this); 
+
+		}
+		/*}}}*/
 };
 
Index: sm/trunk-jpl/src/c/classes/objects/Segment.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/objects/Segment.cpp	(revision 13219)
+++ 	(revision )
@@ -1,86 +1,0 @@
-/*!\file Segment.c
- * \brief: implementation of the Segment object
- */
-
-/*Include files: {{{*/
-#ifdef HAVE_CONFIG_H
-	#include <config.h>
-#else
-#error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!"
-#endif
-
-#include <stdio.h>
-#include <string.h>
-#include "./objects.h"
-#include "../../Container/Container.h"
-#include "../../EnumDefinitions/EnumDefinitions.h"
-#include "../../shared/shared.h"
-#include "../../include/include.h"
-/*}}}*/
-
-/*Segment constructors and destructors:*/
-/*FUNCTION Segment::Segment() default constructor {{{*/
-Segment::Segment(){
-	this->eid=UNDEF;
-	this->x1=UNDEF;
-	this->y1=UNDEF;
-	this->x2=UNDEF;
-	this->y2=UNDEF;
-}
-/*}}}*/
-/*FUNCTION Segment::Segment(int eid, IssmDouble x1,IssmDouble y1,IssmDouble x2, IssmDouble y2){{{*/
-Segment::Segment(int segment_eid, IssmDouble segment_x1,IssmDouble segment_y1,IssmDouble segment_x2, IssmDouble segment_y2){
-
-	this->eid=segment_eid;
-	this->x1=segment_x1;
-	this->y1=segment_y1;
-	this->x2=segment_x2;
-	this->y2=segment_y2;
-
-}
-/*}}}*/
-/*FUNCTION Segment::~Segment(){{{*/
-Segment::~Segment(){
-}
-/*}}}*/
-
-/*Object virtual functions definitions:*/
-/*FUNCTION Segment::Echo{{{*/
-void Segment::Echo(void){
-
-	_printLine_("Segment:");
-	_printLine_("   eid: " << eid);
-	_printLine_("   node 1: " << this->x1 << "|" << this->y1);
-	_printLine_("   node 2: " << this->x2 << "|" << this->y2);
-
-}
-/*}}}*/
-/*FUNCTION Segment::DeepEcho{{{*/
-void Segment::DeepEcho(void){
-	this->Echo();
-}
-/*}}}*/
-/*FUNCTION Segment::Id{{{*/
-int    Segment::Id(void){ return eid; }
-/*}}}*/
-/*FUNCTION Segment::MyRank{{{*/
-int    Segment::MyRank(void){ 
-	extern int my_rank;
-
-	return my_rank; 
-}
-/*}}}*/
-/*FUNCTION Segment::ObjectEnum{{{*/
-int Segment::ObjectEnum(void){
-
-	return SegmentEnum;
-
-}
-/*}}}*/
-/*FUNCTION Segment::copy {{{*/
-Object* Segment::copy() {
-	return new Segment(this->eid,this->x1,this->y1,this->x2,this->y2);
-
-}
-
-/*}}}*/
Index: /issm/trunk-jpl/src/c/classes/objects/Segment.h
===================================================================
--- /issm/trunk-jpl/src/c/classes/objects/Segment.h	(revision 13219)
+++ /issm/trunk-jpl/src/c/classes/objects/Segment.h	(revision 13220)
@@ -11,26 +11,82 @@
 /*}}}*/
 
+template <class doubletype> 
 class Segment: public Object{
 
 	public:
 		int eid;
-		IssmDouble x1;
-		IssmDouble y1;
-		IssmDouble x2;
-		IssmDouble y2;
+		doubletype x1;
+		doubletype y1;
+		doubletype x2;
+		doubletype y2;
 
-		/*Segment constructors, destructors {{{*/
-		Segment();
-		Segment(int eid,IssmDouble x1,IssmDouble y1, IssmDouble x2, IssmDouble y2);
-		~Segment();
+
+		/*Segment constructors, destructors :*/
+		/*FUNCTION Segment() default constructor {{{*/
+		Segment(){
+			this->eid=UNDEF;
+			this->x1=UNDEF;
+			this->y1=UNDEF;
+			this->x2=UNDEF;
+			this->y2=UNDEF;
+		}
 		/*}}}*/
-		/*Object virtual functions definitions:{{{ */
-		void  Echo();
-		void  DeepEcho();
-		int   Id(); 
-		int   MyRank();
-		int   ObjectEnum();
-		Object* copy();
+		/*FUNCTION Segment(int eid, doubletype x1,doubletype y1,doubletype x2, doubletype y2){{{*/
+		Segment(int segment_eid, doubletype segment_x1,doubletype segment_y1,doubletype segment_x2, doubletype segment_y2){
+
+			this->eid=segment_eid;
+			this->x1=segment_x1;
+			this->y1=segment_y1;
+			this->x2=segment_x2;
+			this->y2=segment_y2;
+
+		}
 		/*}}}*/
+		/*FUNCTION ~Segment(){{{*/
+		~Segment(){
+		}
+		/*}}}*/
+
+		/*Object virtual functions definitions:*/
+		/*FUNCTION Echo{{{*/
+		void Echo(void){
+
+			_printLine_("Segment:");
+			_printLine_("   eid: " << eid);
+			_printLine_("   node 1: " << this->x1 << "|" << this->y1);
+			_printLine_("   node 2: " << this->x2 << "|" << this->y2);
+
+		}
+		/*}}}*/
+		/*FUNCTION DeepEcho{{{*/
+		void DeepEcho(void){
+			this->Echo();
+		}
+		/*}}}*/
+		/*FUNCTION Id{{{*/
+		int    Id(void){ return eid; }
+		/*}}}*/
+		/*FUNCTION MyRank{{{*/
+		int    MyRank(void){ 
+			extern int my_rank;
+
+			return my_rank; 
+		}
+		/*}}}*/
+		/*FUNCTION ObjectEnum{{{*/
+		int ObjectEnum(void){
+
+			return SegmentEnum;
+
+		}
+		/*}}}*/
+		/*FUNCTION copy {{{*/
+		Object* copy() {
+			return new Segment(this->eid,this->x1,this->y1,this->x2,this->y2);
+
+		}
+
+		/*}}}*/
+
 };
 
Index: /issm/trunk-jpl/src/c/modules/AddExternalResultx/AddExternalResultx.cpp
===================================================================
--- /issm/trunk-jpl/src/c/modules/AddExternalResultx/AddExternalResultx.cpp	(revision 13219)
+++ /issm/trunk-jpl/src/c/modules/AddExternalResultx/AddExternalResultx.cpp	(revision 13220)
@@ -9,9 +9,9 @@
 #include "../../classes/objects/objects.h"
 		
-void AddExternalResultx( DataSet* results, int enumtype, double value){
+void AddExternalResultx( DataSet* results, int enumtype, IssmDouble value){
 	/* Add new result in into results*/
 	results->AddObject(new DoubleExternalResult(results->Size()+1,enumtype,value,1,0));
 }
-void AddExternalResultx( DataSet* results, int enumtype, double* value, int nraws){
+void AddExternalResultx( DataSet* results, int enumtype, IssmDouble* value, int nraws){
 	/* Add new result in into results*/
 	results->AddObject(new DoubleVecExternalResult(results->Size()+1,enumtype,value,nraws,1,0));
Index: /issm/trunk-jpl/src/c/modules/AddExternalResultx/AddExternalResultx.h
===================================================================
--- /issm/trunk-jpl/src/c/modules/AddExternalResultx/AddExternalResultx.h	(revision 13219)
+++ /issm/trunk-jpl/src/c/modules/AddExternalResultx/AddExternalResultx.h	(revision 13220)
@@ -9,6 +9,6 @@
 
 /* local prototypes: */
-void AddExternalResultx(DataSet* results, int type, double value);
-void AddExternalResultx(DataSet* results, int type, double* value, int nraws);
+void AddExternalResultx(DataSet* results, int type, IssmDouble value);
+void AddExternalResultx(DataSet* results, int type, IssmDouble* value, int nraws);
 
 #endif  /* _ADDEXTERNALRESULT_H */
Index: /issm/trunk-jpl/src/c/modules/ContourToMeshx/ContourToMeshx.cpp
===================================================================
--- /issm/trunk-jpl/src/c/modules/ContourToMeshx/ContourToMeshx.cpp	(revision 13219)
+++ /issm/trunk-jpl/src/c/modules/ContourToMeshx/ContourToMeshx.cpp	(revision 13220)
@@ -11,5 +11,5 @@
 #include "./ContourToMeshx.h"
 
-int ContourToMeshx( Vector<IssmDouble>** pin_nod,Vector<IssmDouble>** pin_elem, double* index, double* x, double* y,DataSet* contours,char* interptype,int nel,int nods, int edgevalue) {
+int ContourToMeshx( Vector<double>** pin_nod,Vector<double>** pin_elem, double* index, double* x, double* y,DataSet* contours,char* interptype,int nel,int nods, int edgevalue) {
 
 	int noerr=1;
@@ -30,9 +30,9 @@
 
 	/*output: */
-	Vector<IssmDouble>* in_nod=NULL;
-	Vector<IssmDouble>* in_elem=NULL;
+	Vector<double>* in_nod=NULL;
+	Vector<double>* in_elem=NULL;
 
-	in_nod=new Vector<IssmDouble>(nods);
-	in_elem=new Vector<IssmDouble>(nel);
+	in_nod=new Vector<double>(nods);
+	in_elem=new Vector<double>(nel);
 
 	/*initialize thread parameters: */
Index: /issm/trunk-jpl/src/c/modules/ContourToMeshx/ContourToMeshx.h
===================================================================
--- /issm/trunk-jpl/src/c/modules/ContourToMeshx/ContourToMeshx.h	(revision 13219)
+++ /issm/trunk-jpl/src/c/modules/ContourToMeshx/ContourToMeshx.h	(revision 13220)
@@ -16,5 +16,5 @@
 	int nods;
 	int edgevalue;
-	Vector<IssmDouble>* in_nod;
+	Vector<double>* in_nod;
 	double* x;
 	double* y;
@@ -24,5 +24,5 @@
 
 /* local prototypes: */
-int ContourToMeshx( Vector<IssmDouble>** pin_nods,Vector<IssmDouble>** pin_elem, double* index, double* x, double* y,DataSet* contours,char* interptype,int nel,int nods, int edgevalue);
+int ContourToMeshx( Vector<double>** pin_nods,Vector<double>** pin_elem, double* index, double* x, double* y,DataSet* contours,char* interptype,int nel,int nods, int edgevalue);
 
 void* ContourToMeshxt(void* vContourToMeshxThreadStruct);
Index: /issm/trunk-jpl/src/c/modules/ContourToMeshx/ContourToMeshxt.cpp
===================================================================
--- /issm/trunk-jpl/src/c/modules/ContourToMeshx/ContourToMeshxt.cpp	(revision 13219)
+++ /issm/trunk-jpl/src/c/modules/ContourToMeshx/ContourToMeshxt.cpp	(revision 13220)
@@ -33,5 +33,5 @@
 	double* x=NULL;
 	double* y=NULL;
-	Vector<IssmDouble>* in_nod=NULL;
+	Vector<double>* in_nod=NULL;
 
 
@@ -55,5 +55,5 @@
 	/*Loop through all contours: */
 	for (i=0;i<contours->Size();i++){
-		Contour* contour=(Contour*)contours->GetObjectByOffset(i);
+		Contour<double>* contour=(Contour<double>*)contours->GetObjectByOffset(i);
 		IsInPoly(in_nod,contour->x,contour->y,contour->nods,x,y,i0,i1,edgevalue);
 	}
Index: /issm/trunk-jpl/src/c/modules/ContourToNodesx/ContourToNodesx.cpp
===================================================================
--- /issm/trunk-jpl/src/c/modules/ContourToNodesx/ContourToNodesx.cpp	(revision 13219)
+++ /issm/trunk-jpl/src/c/modules/ContourToNodesx/ContourToNodesx.cpp	(revision 13220)
@@ -4,5 +4,5 @@
 #include "./ContourToNodesx.h"
 
-int ContourToNodesx( Vector<IssmDouble>** pflags,double* x, double* y, int nods, Contour** contours,int numcontours,int edgevalue){
+int ContourToNodesx( Vector<IssmPDouble>** pflags,double* x, double* y, int nods, Contour<IssmPDouble>** contours,int numcontours,int edgevalue){
 
 	int i;
@@ -10,5 +10,5 @@
 
 	/*Contour:*/
-	Contour* contouri=NULL;
+	Contour<IssmPDouble>* contouri=NULL;
 	int      numnodes;
 	double*  xc=NULL;
@@ -17,7 +17,7 @@
 
 	/*output: */
-	Vector<IssmDouble>* flags=NULL;
+	Vector<IssmPDouble>* flags=NULL;
 
-	flags=new Vector<IssmDouble>(nods);
+	flags=new Vector<IssmPDouble>(nods);
 
 	/*Loop through all contours: */
@@ -39,5 +39,5 @@
 }
 
-int ContourToNodesx( Vector<IssmDouble>** pflags,double* x, double* y, int nods, DataSet* contours, int edgevalue){
+int ContourToNodesx( Vector<IssmPDouble>** pflags,double* x, double* y, int nods, DataSet* contours, int edgevalue){
 
 	int i;
@@ -45,5 +45,5 @@
 
 	/*Contour:*/
-	Contour* contouri=NULL;
+	Contour<IssmPDouble>* contouri=NULL;
 	int      numnodes;
 	double*  xc=NULL;
@@ -52,12 +52,12 @@
 
 	/*output: */
-	Vector<IssmDouble>* flags=NULL;
+	Vector<IssmPDouble>* flags=NULL;
 
-	flags=new Vector<IssmDouble>(nods);
+	flags=new Vector<IssmPDouble>(nods);
 
 	/*Loop through all contours: */
 	if(contours){
 		for (i=0;i<contours->Size();i++){
-			Contour* contour=(Contour*)contours->GetObjectByOffset(i);
+			Contour<IssmPDouble>* contour=(Contour<IssmPDouble>*)contours->GetObjectByOffset(i);
 			IsInPoly(flags,contour->x,contour->y,contour->nods,x,y,0,nods,edgevalue);
 		}
Index: /issm/trunk-jpl/src/c/modules/ContourToNodesx/ContourToNodesx.h
===================================================================
--- /issm/trunk-jpl/src/c/modules/ContourToNodesx/ContourToNodesx.h	(revision 13219)
+++ /issm/trunk-jpl/src/c/modules/ContourToNodesx/ContourToNodesx.h	(revision 13220)
@@ -11,6 +11,6 @@
 
 /* local prototypes: */
-int ContourToNodesx( Vector<IssmDouble>** pflags,double* x, double* y, int nods, Contour** contours,int numcontours,int edgevalue);
-int ContourToNodesx( Vector<IssmDouble>** pflags,double* x, double* y, int nods, DataSet* contours, int edgevalue);
+int ContourToNodesx( Vector<IssmPDouble>** pflags,double* x, double* y, int nods, Contour<IssmPDouble>** contours,int numcontours,int edgevalue);
+int ContourToNodesx( Vector<IssmPDouble>** pflags,double* x, double* y, int nods, DataSet* contours, int edgevalue);
 
 #endif /* _CONTOURTONODESX_H */
Index: /issm/trunk-jpl/src/c/modules/InterpFromGridToMeshx/InterpFromGridToMeshx.cpp
===================================================================
--- /issm/trunk-jpl/src/c/modules/InterpFromGridToMeshx/InterpFromGridToMeshx.cpp	(revision 13219)
+++ /issm/trunk-jpl/src/c/modules/InterpFromGridToMeshx/InterpFromGridToMeshx.cpp	(revision 13220)
@@ -17,8 +17,8 @@
 
 /*InterpFromGridToMeshx{{{*/
-int InterpFromGridToMeshx( Vector<IssmDouble>** pdata_mesh,double* x_in, int x_rows, double* y_in, int y_rows, double* data, int M, int N, double* x_mesh, double* y_mesh, int nods,double default_value, int interpenum){
+int InterpFromGridToMeshx( Vector<IssmPDouble>** pdata_mesh,double* x_in, int x_rows, double* y_in, int y_rows, double* data, int M, int N, double* x_mesh, double* y_mesh, int nods,double default_value, int interpenum){
 
 	/*output: */
-	Vector<IssmDouble>* data_mesh=NULL;
+	Vector<IssmPDouble>* data_mesh=NULL;
 	
 	/*Intermediary*/
@@ -47,5 +47,5 @@
 
 	/*Allocate output vector: */
-	data_mesh=new Vector<IssmDouble>(nods);
+	data_mesh=new Vector<IssmPDouble>(nods);
 
 	/*Find out what kind of coordinates (x_in,y_in) have been given is input*/
@@ -127,5 +127,5 @@
 	double *y             = gate->y;
 	int     nods          = gate->nods;
-	Vector<IssmDouble>*data_mesh     = gate->data_mesh;
+	Vector<IssmPDouble>*data_mesh     = gate->data_mesh;
 	double *data          = gate->data;
 	double  default_value = gate->default_value;
@@ -181,5 +181,5 @@
 					return NULL; /*WARNING: no error because it would blow up the multithreading!*/
 			}
-			if(xIsNan<IssmDouble>(data_value)) data_value=default_value;
+			if(xIsNan<IssmPDouble>(data_value)) data_value=default_value;
 		}
 		else{
Index: /issm/trunk-jpl/src/c/modules/InterpFromGridToMeshx/InterpFromGridToMeshx.h
===================================================================
--- /issm/trunk-jpl/src/c/modules/InterpFromGridToMeshx/InterpFromGridToMeshx.h	(revision 13219)
+++ /issm/trunk-jpl/src/c/modules/InterpFromGridToMeshx/InterpFromGridToMeshx.h	(revision 13220)
@@ -24,8 +24,8 @@
 	double* x_mesh;
 	double* y_mesh;
-	Vector<IssmDouble>*     data_mesh;
+	Vector<IssmPDouble>*     data_mesh;
 } InterpFromGridToMeshxThreadStruct;
 
-int    InterpFromGridToMeshx( Vector<IssmDouble>** pdata_mesh,double* x, int x_rows, double* y, int y_rows, double* data, int M, int N, double* x_mesh, double* y_mesh, int nods, double default_value, int interpenum=BilinearInterpEnum);
+int    InterpFromGridToMeshx( Vector<IssmPDouble>** pdata_mesh,double* x, int x_rows, double* y, int y_rows, double* data, int M, int N, double* x_mesh, double* y_mesh, int nods, double default_value, int interpenum=BilinearInterpEnum);
 void*  InterpFromGridToMeshxt(void* vInterpFromGridToMeshxThreadStruct);
 bool   findindices(int* pn,int* pm,double* x,int x_rows, double* y,int y_rows, double xgrid,double ygrid);
Index: /issm/trunk-jpl/src/c/modules/InterpFromMesh2dx/InterpFromMesh2dx.cpp
===================================================================
--- /issm/trunk-jpl/src/c/modules/InterpFromMesh2dx/InterpFromMesh2dx.cpp	(revision 13219)
+++ /issm/trunk-jpl/src/c/modules/InterpFromMesh2dx/InterpFromMesh2dx.cpp	(revision 13220)
@@ -10,9 +10,9 @@
 #include "../modules.h"
 
-int InterpFromMesh2dx( Vector<IssmDouble>** pdata_prime,double* index_data, double* x_data, double* y_data, int nods_data,int nels_data, double* data, int data_length, double* x_prime, double* y_prime, int nods_prime,
-		double* default_values,int num_default_values,Contour** contours,int numcontours){
+int InterpFromMesh2dx( Vector<IssmPDouble>** pdata_prime,double* index_data, double* x_data, double* y_data, int nods_data,int nels_data, double* data, int data_length, double* x_prime, double* y_prime, int nods_prime,
+		double* default_values,int num_default_values,Contour<IssmPDouble>** contours,int numcontours){
 	
 	/*Output*/
-	Vector<IssmDouble>* data_prime=NULL;
+	Vector<IssmPDouble>* data_prime=NULL;
 
 	/*Intermediary*/
@@ -27,5 +27,5 @@
 
 	/*contours: */
-	Vector<IssmDouble>*    vec_incontour=NULL;
+	Vector<IssmPDouble>*    vec_incontour=NULL;
 	double*    incontour=NULL;
 
@@ -71,5 +71,5 @@
 
 	/*Initialize output*/
-	data_prime=new Vector<IssmDouble>(nods_prime,false,SeqVecType);
+	data_prime=new Vector<IssmPDouble>(nods_prime,false,SeqVecType);
 	if(num_default_values){
 		if(num_default_values==1)for (i=0;i<nods_prime;i++) data_prime->SetValue(i,default_values[0],INS_VAL);
Index: /issm/trunk-jpl/src/c/modules/InterpFromMesh2dx/InterpFromMesh2dx.h
===================================================================
--- /issm/trunk-jpl/src/c/modules/InterpFromMesh2dx/InterpFromMesh2dx.h	(revision 13219)
+++ /issm/trunk-jpl/src/c/modules/InterpFromMesh2dx/InterpFromMesh2dx.h	(revision 13220)
@@ -23,5 +23,5 @@
 	double ymin,ymax;
 	int    nods_prime;
-	Vector<IssmDouble>*    data_prime;
+	Vector<IssmPDouble>*    data_prime;
 	double* x_prime;
 	double* y_prime;
@@ -33,6 +33,6 @@
 } InterpFromMesh2dxThreadStruct;
 
-int InterpFromMesh2dx( Vector<IssmDouble>** pdata_prime,double* index_data, double* x_data, double* y_data, int nods_data,int nels_data, double* data, int data_length, double* x_prime, double* y_prime, int nods_prime,
-		double* default_values,int num_default_values,Contour** contours,int numcontours);
+int InterpFromMesh2dx( Vector<IssmPDouble>** pdata_prime,double* index_data, double* x_data, double* y_data, int nods_data,int nels_data, double* data, int data_length, double* x_prime, double* y_prime, int nods_prime,
+		double* default_values,int num_default_values,Contour<IssmPDouble>** contours,int numcontours);
 
 void* InterpFromMesh2dxt(void* vInterpFromMesh2dxThreadStruct);
Index: /issm/trunk-jpl/src/c/modules/InterpFromMesh2dx/InterpFromMesh2dxt.cpp
===================================================================
--- /issm/trunk-jpl/src/c/modules/InterpFromMesh2dx/InterpFromMesh2dxt.cpp	(revision 13219)
+++ /issm/trunk-jpl/src/c/modules/InterpFromMesh2dx/InterpFromMesh2dxt.cpp	(revision 13220)
@@ -32,5 +32,5 @@
 	double  ymax               = gate->ymax;
 	int     nods_prime         = gate->nods_prime;
-	Vector<IssmDouble>* data_prime         = gate->data_prime;
+	Vector<IssmPDouble>* data_prime         = gate->data_prime;
 	double *x_prime            = gate->x_prime;
 	double *y_prime            = gate->y_prime;
@@ -87,5 +87,5 @@
 						data_value=data[i];
 					}
-					if (xIsNan<IssmDouble>(data_value)){
+					if (xIsNan<IssmPDouble>(data_value)){
 						if(num_default_values==1) data_value=default_values[0];
 						else data_value=default_values[j];
Index: /issm/trunk-jpl/src/c/modules/InterpFromMeshToMesh3dx/InterpFromMeshToMesh3dx.cpp
===================================================================
--- /issm/trunk-jpl/src/c/modules/InterpFromMeshToMesh3dx/InterpFromMeshToMesh3dx.cpp	(revision 13219)
+++ /issm/trunk-jpl/src/c/modules/InterpFromMeshToMesh3dx/InterpFromMeshToMesh3dx.cpp	(revision 13220)
@@ -7,8 +7,8 @@
 #include "../../include/include.h"
 
-int InterpFromMeshToMesh3dx( Vector<IssmDouble>** pdata_prime,double* index_data, double* x_data, double* y_data, double* z_data, int nods_data,int nels_data, double* data, int data_length, double* x_prime, double* y_prime, double* z_prime, int nods_prime,double default_value) {
+int InterpFromMeshToMesh3dx( Vector<IssmPDouble>** pdata_prime,double* index_data, double* x_data, double* y_data, double* z_data, int nods_data,int nels_data, double* data, int data_length, double* x_prime, double* y_prime, double* z_prime, int nods_prime,double default_value) {
 
 	/*Output*/
-	Vector<IssmDouble>* data_prime=NULL;
+	Vector<IssmPDouble>* data_prime=NULL;
 
 	/*Intermediary*/
@@ -54,5 +54,5 @@
 
 	/*Initialize output*/
-	data_prime=new Vector<IssmDouble>(nods_prime);
+	data_prime=new Vector<IssmPDouble>(nods_prime);
 	for (i=0;i<nods_prime;i++) data_prime->SetValue(i,default_value,INS_VAL);
 
@@ -124,5 +124,5 @@
 						data_value=data[i];
 					}
-					if (xIsNan<IssmDouble>(data_value)) data_value=default_value;
+					if (xIsNan<IssmPDouble>(data_value)) data_value=default_value;
 
 					/*insert value and go to the next point*/
Index: /issm/trunk-jpl/src/c/modules/InterpFromMeshToMesh3dx/InterpFromMeshToMesh3dx.h
===================================================================
--- /issm/trunk-jpl/src/c/modules/InterpFromMeshToMesh3dx/InterpFromMeshToMesh3dx.h	(revision 13219)
+++ /issm/trunk-jpl/src/c/modules/InterpFromMeshToMesh3dx/InterpFromMeshToMesh3dx.h	(revision 13220)
@@ -9,5 +9,5 @@
 #include "../../classes/objects/objects.h"
 
-int InterpFromMeshToMesh3dx( Vector<IssmDouble>** pdata_prime,double* index_data, double* x_data, double* y_data, double* z_data, int nods_data,int nels_data, double* data, int data_length, double* x_prime, double* y_prime, double* z_prime, int nods_prime,double default_value);
+int InterpFromMeshToMesh3dx( Vector<IssmPDouble>** pdata_prime,double* index_data, double* x_data, double* y_data, double* z_data, int nods_data,int nels_data, double* data, int data_length, double* x_prime, double* y_prime, double* z_prime, int nods_prime,double default_value);
 
 #endif /* _INTERPFROMMESHTOMESH3DX_H */
Index: sm/trunk-jpl/src/c/modules/MeshPartitionx/MeshPartitionx.cpp
===================================================================
--- /issm/trunk-jpl/src/c/modules/MeshPartitionx/MeshPartitionx.cpp	(revision 13219)
+++ 	(revision )
@@ -1,115 +1,0 @@
-/*!\file:  MeshPartition.cpp
- * \brief partition elements and nodes across a cluster of size numprocs. 
- */
-
-#include "./MeshPartitionx.h"
-#include "../../shared/shared.h"
-#include "../../include/include.h"
-#include "../../toolkits/toolkits.h"
-#include "../../EnumDefinitions/EnumDefinitions.h"
-
-int MeshPartitionx(int** pepart, int** pnpart, int numberofelements,int numberofnodes,IssmDouble* elements,
-		int numberofelements2d,int numberofnodes2d,IssmDouble* elements2d,int numlayers,int elements_width, int dim,int num_procs){
-
-	int noerr=1;
-	int i,j;
-
-	/*Metis partitioning: */
-	int* epart=NULL;
-	int* npart=NULL;
-	int* index=NULL;
-
-	int* epart2d=NULL;
-	int* npart2d=NULL;
-	int* index2d=NULL;
-	int  count=0;
-
-	int  etype=1; //tria mesh see metis/Programs/Io.c
-	int  etype2d=1; //tria mesh see metis/Programs/Io.c
-	int  numflag=0;
-	int  edgecut=1;
-
-	if(dim==2){
-		epart=xNew<int>(numberofelements);
-		npart=xNew<int>(numberofnodes);
-		index=xNew<int>(elements_width*numberofelements);
-		for (i=0;i<numberofelements;i++){
-			for (j=0;j<elements_width;j++){
-				*(index+elements_width*i+j)=reCast<int>(*(elements+elements_width*i+j))-1; //-1 for C indexing in Metis
-			}
-		}
-
-		/*Partition using Metis:*/
-		if (num_procs>1){
-			#ifdef _HAVE_METIS_
-			METIS_PartMeshNodalPatch(&numberofelements,&numberofnodes, index, &etype, &numflag, &num_procs, &edgecut, epart, npart);
-			#endif
-		}
-		else if (num_procs==1){
-			/*METIS does not know how to deal with one cpu only!*/
-			for (i=0;i<numberofelements;i++) epart[i]=0;
-			for (i=0;i<numberofnodes;i++)    npart[i]=0;
-		}
-		else _error_("At least one processor is required");
-	}
-	else{
-		/*We have a 3d mesh, made of a regularly extruded 2d mesh. We first partition the 2d mesh, then we extrude the partition: */
-
-		/*First build concatenated 2d mesh  from 2d_coll and 2d_noncoll: */
-		epart2d=xNew<int>(numberofelements2d);
-		npart2d=xNew<int>(numberofnodes2d); 
-		index2d=xNew<int>(3*numberofelements2d);
-
-		for (i=0;i<numberofelements2d;i++){
-			for (j=0;j<3;j++){
-				*(index2d+3*i+j)=reCast<int>(*(elements2d+3*i+j))-1; //-1 for C indexing in Metis
-			}
-		}
-
-		/*Partition using Metis:*/
-		if (num_procs>1){
-			#ifdef _HAVE_METIS_
-			METIS_PartMeshNodalPatch(&numberofelements2d,&numberofnodes2d, index2d, &etype2d, &numflag, &num_procs, &edgecut, epart2d, npart2d);
-			#endif
-		}
-		else if (num_procs==1){
-			/*METIS does not know how to deal with one cpu only!*/
-			for (i=0;i<numberofelements2d;i++) epart2d[i]=0;
-			for (i=0;i<numberofnodes2d;i++)    npart2d[i]=0;
-		}
-		else _error_("At least one processor is required");
-
-		/*Extrude epart2d to epart, using numlayers: */
-		epart=xNew<int>(numberofelements);
-		
-		count=0;
-		for(i=0;i<(numlayers-1);i++){
-			for(j=0;j<numberofelements2d;j++){
-				epart[count]=epart2d[j];
-				count++;
-			}
-		}
-
-		/*Extrude npart2d to npart, using numlayers: */
-		npart=xNew<int>(numberofnodes);
-		
-		count=0;
-		for(i=0;i<(numlayers);i++){
-			for(j=0;j<numberofnodes2d;j++){
-				npart[count]=npart2d[j];
-				count++;
-			}
-		}
-	}
-	
-	/*Assign output pointer:*/
-	*pepart=epart;
-	*pnpart=npart;
-
-	/*Free ressources: */
-	xDelete<int>(index);
-	xDelete<int>(epart2d);
-	xDelete<int>(npart2d);
-	xDelete<int>(index2d);
-	return noerr;
-}
Index: /issm/trunk-jpl/src/c/modules/MeshPartitionx/MeshPartitionx.h
===================================================================
--- /issm/trunk-jpl/src/c/modules/MeshPartitionx/MeshPartitionx.h	(revision 13219)
+++ /issm/trunk-jpl/src/c/modules/MeshPartitionx/MeshPartitionx.h	(revision 13220)
@@ -7,9 +7,114 @@
 
 #include "../../include/include.h"
+#include "../../EnumDefinitions/EnumDefinitions.h"
+
 
 /* local prototypes: */
-int MeshPartitionx(int** pepart, int** pnpart, int numberofelements,int numberofnodes,IssmDouble* elements,
-		int numberofelements2d,int numberofnodes2d,IssmDouble* elements2d,int numlayers,int elements_width, int dim,int numareas);
+template <class doubletype> 
+int MeshPartitionx(int** pepart, int** pnpart, int numberofelements,int numberofnodes,doubletype* elements,
+		int numberofelements2d,int numberofnodes2d,doubletype* elements2d,int numlayers,int elements_width, int dim,int num_procs){
+
+	int noerr=1;
+	int i,j;
+
+	/*Metis partitioning: */
+	int* epart=NULL;
+	int* npart=NULL;
+	int* index=NULL;
+
+	int* epart2d=NULL;
+	int* npart2d=NULL;
+	int* index2d=NULL;
+	int  count=0;
+
+	int  etype=1; //tria mesh see metis/Programs/Io.c
+	int  etype2d=1; //tria mesh see metis/Programs/Io.c
+	int  numflag=0;
+	int  edgecut=1;
+
+	if(dim==2){
+		epart=xNew<int>(numberofelements);
+		npart=xNew<int>(numberofnodes);
+		index=xNew<int>(elements_width*numberofelements);
+		for (i=0;i<numberofelements;i++){
+			for (j=0;j<elements_width;j++){
+				*(index+elements_width*i+j)=reCast<int>(*(elements+elements_width*i+j))-1; //-1 for C indexing in Metis
+			}
+		}
+
+		/*Partition using Metis:*/
+		if (num_procs>1){
+			#ifdef _HAVE_METIS_
+			METIS_PartMeshNodalPatch(&numberofelements,&numberofnodes, index, &etype, &numflag, &num_procs, &edgecut, epart, npart);
+			#endif
+		}
+		else if (num_procs==1){
+			/*METIS does not know how to deal with one cpu only!*/
+			for (i=0;i<numberofelements;i++) epart[i]=0;
+			for (i=0;i<numberofnodes;i++)    npart[i]=0;
+		}
+		else _error_("At least one processor is required");
+	}
+	else{
+		/*We have a 3d mesh, made of a regularly extruded 2d mesh. We first partition the 2d mesh, then we extrude the partition: */
+
+		/*First build concatenated 2d mesh  from 2d_coll and 2d_noncoll: */
+		epart2d=xNew<int>(numberofelements2d);
+		npart2d=xNew<int>(numberofnodes2d); 
+		index2d=xNew<int>(3*numberofelements2d);
+
+		for (i=0;i<numberofelements2d;i++){
+			for (j=0;j<3;j++){
+				*(index2d+3*i+j)=reCast<int>(*(elements2d+3*i+j))-1; //-1 for C indexing in Metis
+			}
+		}
+
+		/*Partition using Metis:*/
+		if (num_procs>1){
+			#ifdef _HAVE_METIS_
+			METIS_PartMeshNodalPatch(&numberofelements2d,&numberofnodes2d, index2d, &etype2d, &numflag, &num_procs, &edgecut, epart2d, npart2d);
+			#endif
+		}
+		else if (num_procs==1){
+			/*METIS does not know how to deal with one cpu only!*/
+			for (i=0;i<numberofelements2d;i++) epart2d[i]=0;
+			for (i=0;i<numberofnodes2d;i++)    npart2d[i]=0;
+		}
+		else _error_("At least one processor is required");
+
+		/*Extrude epart2d to epart, using numlayers: */
+		epart=xNew<int>(numberofelements);
+		
+		count=0;
+		for(i=0;i<(numlayers-1);i++){
+			for(j=0;j<numberofelements2d;j++){
+				epart[count]=epart2d[j];
+				count++;
+			}
+		}
+
+		/*Extrude npart2d to npart, using numlayers: */
+		npart=xNew<int>(numberofnodes);
+		
+		count=0;
+		for(i=0;i<(numlayers);i++){
+			for(j=0;j<numberofnodes2d;j++){
+				npart[count]=npart2d[j];
+				count++;
+			}
+		}
+	}
 	
+	/*Assign output pointer:*/
+	*pepart=epart;
+	*pnpart=npart;
+
+	/*Free ressources: */
+	xDelete<int>(index);
+	xDelete<int>(epart2d);
+	xDelete<int>(npart2d);
+	xDelete<int>(index2d);
+	return noerr;
+}	
 #endif /* _MESHPARTITIONX_H */
 
Index: /issm/trunk-jpl/src/c/modules/MeshProfileIntersectionx/ElementSegment.cpp
===================================================================
--- /issm/trunk-jpl/src/c/modules/MeshProfileIntersectionx/ElementSegment.cpp	(revision 13219)
+++ /issm/trunk-jpl/src/c/modules/MeshProfileIntersectionx/ElementSegment.cpp	(revision 13220)
@@ -55,5 +55,5 @@
 			yfinal[1]=ysegment[0]+coord2*(ysegment[1]-ysegment[0]);
 
-			segments_dataset->AddObject(new  Segment(el+1,xfinal[0],yfinal[0],xfinal[1],yfinal[1]));
+			segments_dataset->AddObject(new  Segment<double>(el+1,xfinal[0],yfinal[0],xfinal[1],yfinal[1]));
 		}
 		else{
@@ -83,10 +83,10 @@
 		yfinal[1]=ysegment[0]+coord2*(ysegment[1]-ysegment[0]);
 
-		segments_dataset->AddObject(new  Segment(el+1,xfinal[0],yfinal[0],xfinal[1],yfinal[1]));
+		segments_dataset->AddObject(new  Segment<double>(el+1,xfinal[0],yfinal[0],xfinal[1],yfinal[1]));
 	}
 	else{
 		/*No interesections, but the segment might be entirely inside this triangle!: */
 		if ( (NodeInElement(xnodes,ynodes,xsegment[0],ysegment[0])) && (NodeInElement(xnodes,ynodes,xsegment[1],ysegment[1])) ){
-			segments_dataset->AddObject(new  Segment(el+1,xsegment[0],ysegment[0],xsegment[1],ysegment[1]));
+			segments_dataset->AddObject(new  Segment<double>(el+1,xsegment[0],ysegment[0],xsegment[1],ysegment[1]));
 		}
 	}
Index: /issm/trunk-jpl/src/c/modules/MeshProfileIntersectionx/MeshProfileIntersectionx.cpp
===================================================================
--- /issm/trunk-jpl/src/c/modules/MeshProfileIntersectionx/MeshProfileIntersectionx.cpp	(revision 13219)
+++ /issm/trunk-jpl/src/c/modules/MeshProfileIntersectionx/MeshProfileIntersectionx.cpp	(revision 13220)
@@ -4,5 +4,5 @@
 #include "./MeshProfileIntersectionx.h"
 
-void MeshProfileIntersectionx( double** psegments, int* pnumsegs, int* index, double* x, double* y, int nel, int nods,  Contour** contours,int numcontours){
+void MeshProfileIntersectionx( double** psegments, int* pnumsegs, int* index, double* x, double* y, int nel, int nods,  Contour<IssmPDouble>** contours,int numcontours){
 
 	int i,j,k;
@@ -10,5 +10,5 @@
 
 	/*Contour:*/
-	Contour* contouri=NULL;
+	Contour<IssmPDouble>* contouri=NULL;
 	int      numnodes;
 	double*  xc=NULL;
Index: /issm/trunk-jpl/src/c/modules/MeshProfileIntersectionx/MeshProfileIntersectionx.h
===================================================================
--- /issm/trunk-jpl/src/c/modules/MeshProfileIntersectionx/MeshProfileIntersectionx.h	(revision 13219)
+++ /issm/trunk-jpl/src/c/modules/MeshProfileIntersectionx/MeshProfileIntersectionx.h	(revision 13220)
@@ -11,5 +11,5 @@
 
 /* local prototypes: */
-void MeshProfileIntersectionx( double** psegments, int* pnumseg, int* index, double* x, double* y, int nel, int nods,  Contour** contours,int numcontours);
+void MeshProfileIntersectionx( double** psegments, int* pnumseg, int* index, double* x, double* y, int nel, int nods,  Contour<IssmPDouble>** contours,int numcontours);
 void MeshSegmentsIntersection(double** psegments, int* pnumsegs,int* index, double* x, double* y, int nel, int nods, double* xc, double* yc, int numnodes);
 void ElementSegmentsIntersection(DataSet* segments_dataset,int el, double* xnodes,double* ynodes,double* xc,double* yc,int numnodes);
Index: /issm/trunk-jpl/src/c/modules/MeshProfileIntersectionx/MeshSegmentsIntersection.cpp
===================================================================
--- /issm/trunk-jpl/src/c/modules/MeshProfileIntersectionx/MeshSegmentsIntersection.cpp	(revision 13219)
+++ /issm/trunk-jpl/src/c/modules/MeshProfileIntersectionx/MeshSegmentsIntersection.cpp	(revision 13220)
@@ -10,5 +10,5 @@
 	/*output: */
 	double*  segments=NULL;
-	Segment* segment=NULL;
+	Segment<double>* segment=NULL;
 	int     numsegs;
 	
@@ -34,5 +34,5 @@
 	segments=xNew<double>(5*numsegs);
 	for(i=0;i<numsegs;i++){
-		Segment* segment=(Segment*)segments_dataset->GetObjectByOffset(i);
+		Segment<double>* segment=(Segment<double>*)segments_dataset->GetObjectByOffset(i);
 		
 		/*x1,y1,x2,y2 then element_id: */
Index: /issm/trunk-jpl/src/c/modules/PointCloudFindNeighborsx/PointCloudFindNeighborsx.cpp
===================================================================
--- /issm/trunk-jpl/src/c/modules/PointCloudFindNeighborsx/PointCloudFindNeighborsx.cpp	(revision 13219)
+++ /issm/trunk-jpl/src/c/modules/PointCloudFindNeighborsx/PointCloudFindNeighborsx.cpp	(revision 13220)
@@ -4,9 +4,9 @@
 #include "./PointCloudFindNeighborsx.h"
 
-int PointCloudFindNeighborsx( Vector<IssmDouble>** pflags,double* x, double* y, int nods, double mindistance,double multithread){
+int PointCloudFindNeighborsx( Vector<IssmPDouble>** pflags,double* x, double* y, int nods, double mindistance,double multithread){
 
 	/*output: */
-	Vector<IssmDouble>* flags=NULL;
-	flags=new Vector<IssmDouble>(nods);
+	Vector<IssmPDouble>* flags=NULL;
+	flags=new Vector<IssmPDouble>(nods);
 
 	/*threading: */
Index: /issm/trunk-jpl/src/c/modules/PointCloudFindNeighborsx/PointCloudFindNeighborsx.h
===================================================================
--- /issm/trunk-jpl/src/c/modules/PointCloudFindNeighborsx/PointCloudFindNeighborsx.h	(revision 13219)
+++ /issm/trunk-jpl/src/c/modules/PointCloudFindNeighborsx/PointCloudFindNeighborsx.h	(revision 13220)
@@ -11,5 +11,5 @@
 
 /* local prototypes: */
-int PointCloudFindNeighborsx( Vector<IssmDouble>** pflags,double* x, double* y, int nods, double mindistance,double multithread);
+int PointCloudFindNeighborsx( Vector<IssmPDouble>** pflags,double* x, double* y, int nods, double mindistance,double multithread);
 
 /*threading: */
@@ -20,5 +20,5 @@
 	int nods;
 	double mindistance;
-	Vector<IssmDouble>* flags;
+	Vector<IssmPDouble>* flags;
 
 
Index: /issm/trunk-jpl/src/c/modules/PointCloudFindNeighborsx/PointCloudFindNeighborsxt.cpp
===================================================================
--- /issm/trunk-jpl/src/c/modules/PointCloudFindNeighborsx/PointCloudFindNeighborsxt.cpp	(revision 13219)
+++ /issm/trunk-jpl/src/c/modules/PointCloudFindNeighborsx/PointCloudFindNeighborsxt.cpp	(revision 13220)
@@ -17,5 +17,5 @@
 	int     nods;
 	double  mindistance;
-	Vector<IssmDouble>*     flags;
+	Vector<IssmPDouble>*     flags;
 
 	/*recover handle and gate: */
Index: /issm/trunk-jpl/src/c/modules/Reducevectorgtosx/Reducevectorgtosx.cpp
===================================================================
--- /issm/trunk-jpl/src/c/modules/Reducevectorgtosx/Reducevectorgtosx.cpp	(revision 13219)
+++ /issm/trunk-jpl/src/c/modules/Reducevectorgtosx/Reducevectorgtosx.cpp	(revision 13220)
@@ -15,5 +15,5 @@
 	int configuration_type;
 	int ssize;
-	double* yg_serial=NULL;
+	IssmDouble* yg_serial=NULL;
 
 	/*first figure out ssize: */
@@ -51,5 +51,5 @@
 
 	/*Free ressources:*/
-	xDelete<double>(yg_serial);
+	xDelete<IssmDouble>(yg_serial);
 
 	/*Assign output pointers:*/
Index: /issm/trunk-jpl/src/c/modules/TriMeshx/TriMeshx.cpp
===================================================================
--- /issm/trunk-jpl/src/c/modules/TriMeshx/TriMeshx.cpp	(revision 13219)
+++ /issm/trunk-jpl/src/c/modules/TriMeshx/TriMeshx.cpp	(revision 13220)
@@ -20,5 +20,5 @@
 /*}}}*/
 
-void TriMeshx(Matrix<IssmDouble>** pindex,Vector<IssmDouble>** px,Vector<IssmDouble>** py,Matrix<IssmDouble>** psegments,Vector<IssmDouble>** psegmentmarkerlist,DataSet* domain,DataSet* rifts,double area){
+void TriMeshx(Matrix<IssmPDouble>** pindex,Vector<IssmPDouble>** px,Vector<IssmPDouble>** py,Matrix<IssmPDouble>** psegments,Vector<IssmPDouble>** psegmentmarkerlist,DataSet* domain,DataSet* rifts,double area){
 
 	/*indexing: */
@@ -26,15 +26,15 @@
 
 	/*output: */
-	double* index=NULL;
-	Matrix<IssmDouble>* index_matrix=NULL;
+	IssmPDouble* index=NULL;
+	Matrix<IssmPDouble>* index_matrix=NULL;
 	double* x=NULL;
 	double* y=NULL;
 	double* segments=NULL;
-	Matrix<IssmDouble>* segments_matrix=NULL;
+	Matrix<IssmPDouble>* segments_matrix=NULL;
 	double* segmentmarkerlist=NULL;
 
 	/*intermediary: */
 	int      counter,counter2,backcounter;
-	Contour* contour=NULL;
+	Contour<IssmPDouble>* contour=NULL;
 
 	/* Triangle structures needed to call Triangle library routines: */
@@ -45,9 +45,9 @@
 	in.numberofpoints=0;
 	for (i=0;i<domain->Size();i++){
-		contour=(Contour*)domain->GetObjectByOffset(i);
+		contour=(Contour<IssmPDouble>*)domain->GetObjectByOffset(i);
 		in.numberofpoints+=contour->nods-1;
 	}
 	for (i=0;i<rifts->Size();i++){
-		contour=(Contour*)rifts->GetObjectByOffset(i);
+		contour=(Contour<IssmPDouble>*)rifts->GetObjectByOffset(i);
 		in.numberofpoints+=contour->nods;
 	}
@@ -61,5 +61,5 @@
 	counter=0;
 	for (i=0;i<domain->Size();i++){
-		contour=(Contour*)domain->GetObjectByOffset(i);
+		contour=(Contour<IssmPDouble>*)domain->GetObjectByOffset(i);
 		for (j=0;j<contour->nods-1;j++){
 			in.pointlist[2*counter+0]=contour->x[j];
@@ -69,5 +69,5 @@
 	}
 	for (i=0;i<rifts->Size();i++){
-		contour=(Contour*)rifts->GetObjectByOffset(i);
+		contour=(Contour<IssmPDouble>*)rifts->GetObjectByOffset(i);
 		for (j=0;j<contour->nods;j++){
 			in.pointlist[2*counter+0]=contour->x[j];
@@ -88,9 +88,9 @@
 	in.numberofsegments=0;
 	for (i=0;i<domain->Size();i++){
-		contour=(Contour*)domain->GetObjectByOffset(i);
+		contour=(Contour<IssmPDouble>*)domain->GetObjectByOffset(i);
 		in.numberofsegments+=contour->nods-1;
 	}
 	for(i=0;i<rifts->Size();i++){
-		contour=(Contour*)rifts->GetObjectByOffset(i);
+		contour=(Contour<IssmPDouble>*)rifts->GetObjectByOffset(i);
 		/*for rifts, we have one less segment as we have vertices*/
 		in.numberofsegments+=contour->nods-1;
@@ -102,5 +102,5 @@
 	backcounter=0;
 	for (i=0;i<domain->Size();i++){
-		contour=(Contour*)domain->GetObjectByOffset(i);
+		contour=(Contour<IssmPDouble>*)domain->GetObjectByOffset(i);
 		for (j=0;j<contour->nods-2;j++){
 			in.segmentlist[2*counter+0]=counter;
@@ -118,5 +118,5 @@
 	counter2=counter;
 	for (i=0;i<rifts->Size();i++){
-		contour=(Contour*)rifts->GetObjectByOffset(i);
+		contour=(Contour<IssmPDouble>*)rifts->GetObjectByOffset(i);
 		for (j=0;j<(contour->nods-1);j++){
 			in.segmentlist[2*counter2+0]=counter;
@@ -137,5 +137,5 @@
 		in.holelist = xNew<REAL>(in.numberofholes*2);
 		for (i=0;i<domain->Size()-1;i++){
-			contour=(Contour*)domain->GetObjectByOffset(i+1);
+			contour=(Contour<IssmPDouble>*)domain->GetObjectByOffset(i+1);
 			GridInsideHole(&in.holelist[2*i+0],&in.holelist[2*i+1],contour->nods-1,contour->x,contour->y);
 		}
@@ -192,12 +192,12 @@
 
 	/*Output : */
-	index_matrix=new Matrix<IssmDouble>(index,out.numberoftriangles,3,1,SeqMatType);
+	index_matrix=new Matrix<IssmPDouble>(index,out.numberoftriangles,3,1.0,SeqMatType);
 	*pindex=index_matrix;
 	
-	segments_matrix=new Matrix<IssmDouble>(segments,out.numberofsegments,3,1,SeqMatType);
+	segments_matrix=new Matrix<IssmPDouble>(segments,out.numberofsegments,3,1.0,SeqMatType);
 	*psegments=segments_matrix;
 
-	*px=new Vector<IssmDouble>(x,out.numberofpoints,SeqMatType);
-	*py=new Vector<IssmDouble>(y,out.numberofpoints,SeqMatType);
-	*psegmentmarkerlist=new Vector<IssmDouble>(segmentmarkerlist,out.numberofsegments,SeqMatType);
+	*px=new Vector<IssmPDouble>(x,out.numberofpoints,SeqMatType);
+	*py=new Vector<IssmPDouble>(y,out.numberofpoints,SeqMatType);
+	*psegmentmarkerlist=new Vector<IssmPDouble>(segmentmarkerlist,out.numberofsegments,SeqMatType);
 }
Index: /issm/trunk-jpl/src/c/modules/TriMeshx/TriMeshx.h
===================================================================
--- /issm/trunk-jpl/src/c/modules/TriMeshx/TriMeshx.h	(revision 13219)
+++ /issm/trunk-jpl/src/c/modules/TriMeshx/TriMeshx.h	(revision 13220)
@@ -11,5 +11,5 @@
 
 /* local prototypes: */
-void TriMeshx(Matrix<IssmDouble>** pindex,Vector<IssmDouble>** px,Vector<IssmDouble>** py,Matrix<IssmDouble>** psegments,Vector<IssmDouble>** psegmentmarkerlist,DataSet* domain,DataSet* rifts,double area);
+void TriMeshx(Matrix<IssmPDouble>** pindex,Vector<IssmPDouble>** px,Vector<IssmPDouble>** py,Matrix<IssmPDouble>** psegments,Vector<IssmPDouble>** psegmentmarkerlist,DataSet* domain,DataSet* rifts,double area);
 
 #endif  /* _TRIMESHX_H */
Index: sm/trunk-jpl/src/c/shared/Exp/DomainOutlineRead.cpp
===================================================================
--- /issm/trunk-jpl/src/c/shared/Exp/DomainOutlineRead.cpp	(revision 13219)
+++ 	(revision )
@@ -1,153 +1,0 @@
-/*!\file:  DomainOutlineRead.cpp
- * \brief DomainOutlineRead.c: read the vertex coordinates defined in a domain 
- * outline from Argus (.exp file). The first contour in the file is for 
- * the outside domain outline. The following contours represent holes in
- * the domain.
- */
-
-#include <stdio.h>
-#include <cstring> 
-#include "../Alloc/alloc.h"
-#include "../../include/include.h"
-#include "../../classes/classes.h"
-#include "../Exceptions/exceptions.h"
-#include "../../Container/DataSet.h"
-
-int DomainOutlineRead(int* pnprof,int** pprofnvertices,double*** ppprofx,double*** ppprofy,bool** pclosed,char* domainname){
-	
-	/*indexing: */
-	int i,counter;
-
-	/*I/O: */
-	FILE   *fid = NULL;
-	char    chardummy[256];
-	double  ddummy;
-
-	/*output: */
-	int      nprof;                //number of profiles in the domainname file
-	int     *profnvertices = NULL; //array holding the number of vertices for the nprof profiles
-	double **pprofx        = NULL; //array of profiles x coordinates
-	double **pprofy        = NULL; //array of profiles y coordinates
-	bool    *closed        = NULL; //array holding closed flags for the nprof profiles
-
-	/*For each profile: */
-	int     n;
-	double *x  = NULL;
-	double *y  = NULL;
-	bool    cl;
-
-	/*open domain outline file for reading: */
-	if ((fid=fopen(domainname,"r"))==NULL){
-		_error_("could not find domain file " << domainname); 
-	}
-
-	/*Do a first pass through the domainname file, to figure out how many profiles 
-	 *we need to read: */
-	nprof=1;
-	for(;;){
-		fscanf(fid,"%256s %256s\n",chardummy,chardummy);
-		fscanf(fid,"%256s %256s\n",chardummy,chardummy);
-		fscanf(fid,"%256s %256s %256s %256s\n",chardummy,chardummy,chardummy,chardummy);
-		fscanf(fid,"%20u %256s\n",&n,chardummy);
-		fscanf(fid,"%256s %256s %256s %256s %256s\n",chardummy,chardummy,chardummy,chardummy,chardummy);
-		for (i=0;i<n;i++){
-			fscanf(fid,"%20lf %20lf\n",&ddummy,&ddummy);
-		}
-		/*Ok, we have faked one profile reading, check whether we are at the end of the file, otherwise, keep fake reading next profile:*/
-		if (feof(fid)){
-			break;
-		}
-		nprof++;
-	}
-	
-	/*Allocate and initialize all the profiles: */
-	profnvertices=xNew<int>(nprof);
-	pprofx=xNew<double*>(nprof);
-	pprofy=xNew<double*>(nprof);
-	for (i=0;i<nprof;i++){
-		pprofx[i]=NULL;
-		pprofy[i]=NULL;
-	}
-	closed=xNew<bool>(nprof);
-
-	/*Reaset file pointer to beginning of file: */
-	fseek(fid,0,SEEK_SET);
-
-	/*Start reading profiles: */
-	for(counter=0;counter<nprof;counter++){
-
-		/*Skip header: */
-		fscanf(fid,"%256s %256s\n",chardummy,chardummy);
-		fscanf(fid,"%256s %256s\n",chardummy,chardummy);
-		fscanf(fid,"%256s %256s %256s %256s\n",chardummy,chardummy,chardummy,chardummy);
-		
-		/*Get number of profile vertices: */
-		fscanf(fid,"%20u %256s\n",&n,chardummy);
-	
-		/*Skip next line: */
-		fscanf(fid,"%256s %256s %256s %256s %256s\n",chardummy,chardummy,chardummy,chardummy,chardummy);
-
-		/*Allocate vertices: */
-		x=xNew<double>(n);
-		y=xNew<double>(n);
-		
-		/*Read vertices: */
-		for (i=0;i<n;i++){
-			fscanf(fid,"%20lf %20lf\n",&x[i],&y[i]);
-		}
-
-		/*Now check that we are dealing with open contours: */
-		cl=false;
-		if((x[0]==x[n-1]) && (y[0]==y[n-1])){
-			cl=true;
-		}
-
-		/*Assign pointers: */
-		profnvertices[counter]=n;
-		pprofx[counter]=x;
-		pprofy[counter]=y;
-		closed[counter]=cl;
-	}
-
-	/*close domain outline file: */
-	fclose(fid);
-
-	/*Assign output pointers: */
-	*pnprof=nprof;
-	*pprofnvertices=profnvertices;
-	*ppprofx=pprofx;
-	*ppprofy=pprofy;
-	if(pclosed)
-	 *pclosed=closed;
-	else
-	 xDelete<bool>(closed);
-}
-
-DataSet* DomainOutlineRead(char* domainname){
-
-	/*intermediary: */
-	int       nprof;
-	int      *profnvertices = NULL;
-	double  **pprofx        = NULL;
-	double  **pprofy        = NULL;
-	Contour  *contour       = NULL;
-
-	/*output: */
-	DataSet* domain=NULL;
-
-	/*If domainname is an empty string, return empty dataset*/
-	if (strcmp(domainname,"")==0){
-		nprof=0;
-	}
-	else{
-		DomainOutlineRead(&nprof,&profnvertices,&pprofx, &pprofy, NULL,domainname);
-	}
-
-	/*now create dataset of contours: */
-	domain=new DataSet(0);
-
-	for(int i=0;i<nprof;i++){
-		domain->AddObject(new Contour(i,profnvertices[i],pprofx[i],pprofy[i],1));
-	}
-	return domain;
-}
Index: /issm/trunk-jpl/src/c/shared/Exp/IsInPoly.cpp
===================================================================
--- /issm/trunk-jpl/src/c/shared/Exp/IsInPoly.cpp	(revision 13219)
+++ /issm/trunk-jpl/src/c/shared/Exp/IsInPoly.cpp	(revision 13220)
@@ -14,86 +14,4 @@
 #endif
 
-/*IsInPoly {{{*/
-int IsInPoly(Vector<IssmDouble>* in,double* xc,double* yc,int numvertices,double* x,double* y,int i0,int i1, int edgevalue){
-
-	int i;
-	double x0,y0;
-	double value;
-	double xmin=xc[0];
-	double xmax=xc[0];
-	double ymin=yc[0];
-	double ymax=yc[0];
-
-	/*Get extrema*/
-	for (i=1;i<numvertices;i++){
-		if(xc[i]<xmin) xmin=xc[i];
-		if(xc[i]>xmax) xmax=xc[i];
-		if(yc[i]<ymin) ymin=yc[i];
-		if(yc[i]>ymax) ymax=yc[i];
-	}
-
-	/*Go through all vertices of the mesh:*/
-	for (i=i0;i<i1;i++){
-
-		//Get current value of value[i] -> do not change it if != 0
-		in->GetValue(&value,i);
-		if (value){
-			/*this vertex already is inside one of the contours, continue*/
-			continue;
-		}
-
-		/*pick up vertex (x[i],y[i]) and figure out if located inside contour (xc,yc)*/
-		x0=x[i]; y0=y[i];
-		if(x0<xmin || x0>xmax || y0<ymin || y0>ymax){
-			value=0;
-		}
-		else{
-			value=pnpoly(numvertices,xc,yc,x0,y0,edgevalue);
-		}
-		in->SetValue(i,value,INS_VAL);
-	}
-	 return 1;
-}/*}}}*/
-/*IsOutsidePoly {{{*/
-int IsOutsidePoly(Vector<IssmDouble>* in,double* xc,double* yc,int numvertices,double* x,double* y,int i0,int i1, int edgevalue){
-
-	int i,j;
-	double x0,y0;
-	double value;
-	double xmin=xc[0];
-	double xmax=xc[0];
-	double ymin=yc[0];
-	double ymax=yc[0];
-
-	/*Get extrema*/
-	for (i=1;i<numvertices;i++){
-		if(xc[i]<xmin) xmin=xc[i];
-		if(xc[i]>xmax) xmax=xc[i];
-		if(yc[i]<ymin) ymin=yc[i];
-		if(yc[i]>ymax) ymax=yc[i];
-	}
-
-	/*Go through all vertices of the mesh:*/
-	for (i=i0;i<i1;i++){
-
-		//Get current value of value[i] -> do not change it if != 0
-		in->GetValue(&value,i);
-		if (value){
-			/*this vertex already is inside one of the contours, continue*/
-			continue;
-		}
-
-		/*pick up vertex (x[i],y[i]) and figure out if located inside contour (xc,yc)*/
-		x0=x[i]; y0=y[i];
-		if(x0<xmin || x0>xmax || y0<ymin || y0>ymax){
-			value=1;
-		}
-		else{
-			value=1-pnpoly(numvertices,xc,yc,x0,y0,edgevalue);
-		}
-		in->SetValue(i,value,INS_VAL);
-	}
-	return 1;
-}/*}}}*/
 /*pnpoly{{{*/
 int pnpoly(int npol, double *xp, double *yp, double x, double y, int edgevalue) {
Index: /issm/trunk-jpl/src/c/shared/Exp/exp.h
===================================================================
--- /issm/trunk-jpl/src/c/shared/Exp/exp.h	(revision 13219)
+++ /issm/trunk-jpl/src/c/shared/Exp/exp.h	(revision 13220)
@@ -7,15 +7,242 @@
 #define _EXP_H_
 
-#include "../../classes/objects/objects.h"
-#include "../../toolkits/toolkits.h"
-
-int IsInPoly(Vector<IssmDouble>* in,double* xc,double* yc,int numvertices,double* x,double* y,int i0,int i1, int edgevalue);
-int IsOutsidePoly(Vector<IssmDouble>* in,double* xc,double* yc,int numvertices,double* x,double* y,int i0,int i1, int edgevalue);
+#include "../../classes/objects/Contour.h"
+#include "../../shared/Numerics/recast.h"
+#include "../../Container/Container.h"
+
 int IsInPolySerial(double* in,double* xc,double* yc,int numvertices,double* x,double* y,int nods, int edgevalue);
 int DomainOutlineWrite(int nprof,int* profnvertices,double** pprofx,double** pprofy,bool* closed,char* domainname);
 int pnpoly(int npol, double *xp, double *yp, double x, double y, int edgevalue);
 
-int      DomainOutlineRead(int* pnprof,int** pprofnvertices,double*** ppprofx,double*** ppprofy,bool** pclosed,char* domainname);
-DataSet* DomainOutlineRead(char* domainname);
+
+/*IsInPoly {{{*/
+template <class doubletype>
+int IsInPoly(Vector<doubletype>* in,double* xc,double* yc,int numvertices,double* x,double* y,int i0,int i1, int edgevalue){
+
+	int i;
+	double x0,y0;
+	doubletype value;
+	double xmin=xc[0];
+	double xmax=xc[0];
+	double ymin=yc[0];
+	double ymax=yc[0];
+
+	/*Get extrema*/
+	for (i=1;i<numvertices;i++){
+		if(xc[i]<xmin) xmin=xc[i];
+		if(xc[i]>xmax) xmax=xc[i];
+		if(yc[i]<ymin) ymin=yc[i];
+		if(yc[i]>ymax) ymax=yc[i];
+	}
+
+	/*Go through all vertices of the mesh:*/
+	for (i=i0;i<i1;i++){
+
+		//Get current value of value[i] -> do not change it if != 0
+		in->GetValue(&value,i);
+		if (reCast<bool,doubletype>(value)){
+			/*this vertex already is inside one of the contours, continue*/
+			continue;
+		}
+
+		/*pick up vertex (x[i],y[i]) and figure out if located inside contour (xc,yc)*/
+		x0=x[i]; y0=y[i];
+		if(x0<xmin || x0>xmax || y0<ymin || y0>ymax){
+			value=0;
+		}
+		else{
+			value=pnpoly(numvertices,xc,yc,x0,y0,edgevalue);
+		}
+		in->SetValue(i,value,INS_VAL);
+	}
+	 return 1;
+}/*}}}*/
+/*IsOutsidePoly {{{*/
+template <class doubletype>
+int IsOutsidePoly(Vector<doubletype>* in,double* xc,double* yc,int numvertices,double* x,double* y,int i0,int i1, int edgevalue){
+
+	int i,j;
+	double x0,y0;
+	doubletype value;
+	double xmin=xc[0];
+	double xmax=xc[0];
+	double ymin=yc[0];
+	double ymax=yc[0];
+
+	/*Get extrema*/
+	for (i=1;i<numvertices;i++){
+		if(xc[i]<xmin) xmin=xc[i];
+		if(xc[i]>xmax) xmax=xc[i];
+		if(yc[i]<ymin) ymin=yc[i];
+		if(yc[i]>ymax) ymax=yc[i];
+	}
+
+	/*Go through all vertices of the mesh:*/
+	for (i=i0;i<i1;i++){
+
+		//Get current value of value[i] -> do not change it if != 0
+		in->GetValue(&value,i);
+		if (reCast<bool,doubletype>(value)){
+			/*this vertex already is inside one of the contours, continue*/
+			continue;
+		}
+
+		/*pick up vertex (x[i],y[i]) and figure out if located inside contour (xc,yc)*/
+		x0=x[i]; y0=y[i];
+		if(x0<xmin || x0>xmax || y0<ymin || y0>ymax){
+			value=1;
+		}
+		else{
+			value=1-pnpoly(numvertices,xc,yc,x0,y0,edgevalue);
+		}
+		in->SetValue(i,value,INS_VAL);
+	}
+	return 1;
+}/*}}}*/
+/*DomainOutlineRead{{{*/
+template <class doubletype>
+int DomainOutlineRead(int* pnprof,int** pprofnvertices,doubletype*** ppprofx,doubletype*** ppprofy,bool** pclosed,char* domainname){
+	
+	/*indexing: */
+	int i,counter;
+
+	/*I/O: */
+	FILE   *fid = NULL;
+	char    chardummy[256];
+	double  ddummy;
+
+	/*output: */
+	int      nprof;                //number of profiles in the domainname file
+	int     *profnvertices = NULL; //array holding the number of vertices for the nprof profiles
+	doubletype **pprofx        = NULL; //array of profiles x coordinates
+	doubletype **pprofy        = NULL; //array of profiles y coordinates
+	bool    *closed        = NULL; //array holding closed flags for the nprof profiles
+
+	/*For each profile: */
+	int     n;
+	doubletype *x  = NULL;
+	doubletype *y  = NULL;
+	double  xi,yi;
+	bool    cl;
+
+	/*open domain outline file for reading: */
+	if ((fid=fopen(domainname,"r"))==NULL){
+		_error_("could not find domain file " << domainname); 
+	}
+
+	/*Do a first pass through the domainname file, to figure out how many profiles 
+	 *we need to read: */
+	nprof=1;
+	for(;;){
+		fscanf(fid,"%256s %256s\n",chardummy,chardummy);
+		fscanf(fid,"%256s %256s\n",chardummy,chardummy);
+		fscanf(fid,"%256s %256s %256s %256s\n",chardummy,chardummy,chardummy,chardummy);
+		fscanf(fid,"%20u %256s\n",&n,chardummy);
+		fscanf(fid,"%256s %256s %256s %256s %256s\n",chardummy,chardummy,chardummy,chardummy,chardummy);
+		for (i=0;i<n;i++){
+			fscanf(fid,"%20lf %20lf\n",&ddummy,&ddummy);
+		}
+		/*Ok, we have faked one profile reading, check whether we are at the end of the file, otherwise, keep fake reading next profile:*/
+		if (feof(fid)){
+			break;
+		}
+		nprof++;
+	}
+	
+	/*Allocate and initialize all the profiles: */
+	profnvertices=xNew<int>(nprof);
+	pprofx=xNew<doubletype*>(nprof);
+	pprofy=xNew<doubletype*>(nprof);
+	for (i=0;i<nprof;i++){
+		pprofx[i]=NULL;
+		pprofy[i]=NULL;
+	}
+	closed=xNew<bool>(nprof);
+
+	/*Reaset file pointer to beginning of file: */
+	fseek(fid,0,SEEK_SET);
+
+	/*Start reading profiles: */
+	for(counter=0;counter<nprof;counter++){
+
+		/*Skip header: */
+		fscanf(fid,"%256s %256s\n",chardummy,chardummy);
+		fscanf(fid,"%256s %256s\n",chardummy,chardummy);
+		fscanf(fid,"%256s %256s %256s %256s\n",chardummy,chardummy,chardummy,chardummy);
+		
+		/*Get number of profile vertices: */
+		fscanf(fid,"%20u %256s\n",&n,chardummy);
+	
+		/*Skip next line: */
+		fscanf(fid,"%256s %256s %256s %256s %256s\n",chardummy,chardummy,chardummy,chardummy,chardummy);
+
+		/*Allocate vertices: */
+		x=xNew<doubletype>(n);
+		y=xNew<doubletype>(n);
+		
+		/*Read vertices: */
+		for (i=0;i<n;i++){
+			fscanf(fid,"%20lf %20lf\n",&xi,&yi);
+			x[i]=xi; y[i]=yi;
+		}
+
+		/*Now check that we are dealing with open contours: */
+		cl=false;
+		if((x[0]==x[n-1]) && (y[0]==y[n-1])){
+			cl=true;
+		}
+
+		/*Assign pointers: */
+		profnvertices[counter]=n;
+		pprofx[counter]=x;
+		pprofy[counter]=y;
+		closed[counter]=cl;
+	}
+
+	/*close domain outline file: */
+	fclose(fid);
+
+	/*Assign output pointers: */
+	*pnprof=nprof;
+	*pprofnvertices=profnvertices;
+	*ppprofx=pprofx;
+	*ppprofy=pprofy;
+	if(pclosed)
+	 *pclosed=closed;
+	else
+	 xDelete<bool>(closed);
+} /*}}}*/
+/*DataSet* DomainOutlineRead(char* domainname){{{*/
+template <class doubletype>
+DataSet* DomainOutlineRead(char* domainname){
+
+	/*intermediary: */
+	int       nprof;
+	int      *profnvertices = NULL;
+	doubletype  **pprofx        = NULL;
+	doubletype  **pprofy        = NULL;
+	Contour<doubletype>* contour       = NULL;
+
+	/*output: */
+	DataSet* domain=NULL;
+
+	/*If domainname is an empty string, return empty dataset*/
+	if (strcmp(domainname,"")==0){
+		nprof=0;
+	}
+	else{
+		DomainOutlineRead(&nprof,&profnvertices,&pprofx, &pprofy, NULL,domainname);
+	}
+
+	/*now create dataset of contours: */
+	domain=new DataSet(0);
+
+	for(int i=0;i<nprof;i++){
+		domain->AddObject(new Contour<doubletype>(i,profnvertices[i],pprofx[i],pprofy[i],1));
+	}
+	return domain;
+} /*}}}*/
+
+
 
 
Index: /issm/trunk-jpl/src/modules/ContourToMesh/ContourToMesh.cpp
===================================================================
--- /issm/trunk-jpl/src/modules/ContourToMesh/ContourToMesh.cpp	(revision 13219)
+++ /issm/trunk-jpl/src/modules/ContourToMesh/ContourToMesh.cpp	(revision 13220)
@@ -63,5 +63,5 @@
 	FetchData(&edgevalue,EDGEVALUE);
 	FetchData(&contourname,CONTOURNAME);
-	contours=DomainOutlineRead(contourname);
+	contours=DomainOutlineRead<double>(contourname);
 
 	/*Fetch  interptype: */
Index: /issm/trunk-jpl/src/modules/ContourToNodes/ContourToNodes.cpp
===================================================================
--- /issm/trunk-jpl/src/modules/ContourToNodes/ContourToNodes.cpp	(revision 13219)
+++ /issm/trunk-jpl/src/modules/ContourToNodes/ContourToNodes.cpp	(revision 13220)
@@ -41,6 +41,6 @@
 	mxArray*  matlabstructure=NULL;
 	int numcontours;
-	Contour** contours=NULL;
-	Contour*  contouri=NULL;
+	Contour<double>** contours=NULL;
+	Contour<double>*  contouri=NULL;
 
 	/*Boot module: */
@@ -68,8 +68,8 @@
 
 	numcontours=mxGetNumberOfElements(matlabstructure);
-	contours=xNew<Contour*>(numcontours);
+	contours=xNew<Contour<double>*>(numcontours);
 	for(i=0;i<numcontours;i++){
 		//allocate
-		contouri=xNew<Contour>(1);
+		contouri=xNew<Contour<double> >(1);
 		//retrieve dimension of this contour.
 		contouri->nods=(int)mxGetScalar(mxGetField(matlabstructure,i,"nods"));
Index: /issm/trunk-jpl/src/modules/InterpFromMesh2d/InterpFromMesh2d.cpp
===================================================================
--- /issm/trunk-jpl/src/modules/InterpFromMesh2d/InterpFromMesh2d.cpp	(revision 13219)
+++ /issm/trunk-jpl/src/modules/InterpFromMesh2d/InterpFromMesh2d.cpp	(revision 13220)
@@ -50,7 +50,7 @@
 	//contours
 	mxArray*  matlabstructure=NULL;
-	Contour** contours=NULL;
+	Contour<double>** contours=NULL;
 	int       numcontours;
-	Contour*  contouri=NULL;
+	Contour<double>*  contouri=NULL;
 	int       i;
 
@@ -100,8 +100,8 @@
 		/*contours: */
 		numcontours=mxGetNumberOfElements(matlabstructure);
-		contours=xNew<Contour*>(numcontours);
+		contours=xNew<Contour<double> *>(numcontours);
 		for(i=0;i<numcontours;i++){
 			//allocate
-			contouri=xNew<Contour>(1);
+			contouri=xNew<Contour<double> >(1);
 			//retrieve dimension of this contour.
 			contouri->nods=(int)mxGetScalar(mxGetField(matlabstructure,i,"nods"));
Index: /issm/trunk-jpl/src/modules/MeshProfileIntersection/MeshProfileIntersection.cpp
===================================================================
--- /issm/trunk-jpl/src/modules/MeshProfileIntersection/MeshProfileIntersection.cpp	(revision 13219)
+++ /issm/trunk-jpl/src/modules/MeshProfileIntersection/MeshProfileIntersection.cpp	(revision 13220)
@@ -34,7 +34,7 @@
 	//contours
 	mxArray*  matlabstructure=NULL;
-	Contour** contours=NULL;
+	Contour<double>** contours=NULL;
 	int       numcontours;
-	Contour*  contouri=NULL;
+	Contour<double>*  contouri=NULL;
 
 	/* output: */
@@ -67,8 +67,8 @@
 	//contours
 	numcontours=mxGetNumberOfElements(matlabstructure);
-	contours=xNew<Contour*>(numcontours);
+	contours=xNew<Contour<double>*>(numcontours);
 	for(i=0;i<numcontours;i++){
 		//allocate
-		contouri=xNew<Contour>(1);
+		contouri=xNew<Contour<double> >(1);
 		//retrieve dimension of this contour.
 		contouri->nods=(int)mxGetScalar(mxGetField(matlabstructure,i,"nods"));
Index: /issm/trunk-jpl/src/modules/TriMesh/TriMesh.cpp
===================================================================
--- /issm/trunk-jpl/src/modules/TriMesh/TriMesh.cpp	(revision 13219)
+++ /issm/trunk-jpl/src/modules/TriMesh/TriMesh.cpp	(revision 13220)
@@ -35,6 +35,6 @@
 
 	/*Read domain outline: */
-	domain = DomainOutlineRead(domainname);
-	rifts  = DomainOutlineRead(riftsname);
+	domain = DomainOutlineRead<double>(domainname);
+	rifts  = DomainOutlineRead<double>(riftsname);
 
 	/*call x core: */