diff --git a/tin.h b/tin.h
index f2e81c4..e563bd3 100644
--- a/tin.h
+++ b/tin.h
@@ -13,6 +13,8 @@
 #include "vector"
 #include "algorithm"
 
+#include "iostream"
+
 #define ZERO 1e-5
 
 // Start vertex definition
@@ -77,11 +79,14 @@ bool operator ==(const edge &a, const edge &b) // overload the == operator for e
 // End edge definition
 
 // Start triangle definition
+struct dem_point;
+
 struct triangle
 {
 	vertex2dc *vert[3]; // vertex of the triangle
 	double cx, cy; // center of the triangle's circumcircle
 	double cr; // radius of the circumcircle
+	std::vector<dem_point*> circum_dem;
 
 	triangle() {vert[0] = vert[1] = vert[2] = nullptr;}
 	triangle(vertex2dc *v0ptr, vertex2dc *v1ptr, vertex2dc *v2ptr) {set(v0ptr, v1ptr, v2ptr);}
@@ -264,6 +269,7 @@ void dem2tin(const std::vector<double> &dem, double xmin, double xmax, double ym
 			if ((dist - out_tris[t]->cr) <= ZERO) // Points on the circumcircle are also included
 			{
 				dem_grid[i]->circum_host.push_back(out_tris[t]);
+				out_tris[t]->circum_dem.push_back(dem_grid[i]);
 				// no beak here. There might be more than one triangle's circumcircle includes the DEM location
 			}
 		}
@@ -301,12 +307,12 @@ void dem2tin(const std::vector<double> &dem, double xmin, double xmax, double ym
 			cnst_tri.push_back(dem_grid[0]->circum_host[i]);
 		}
 
-		for (int i = 0; i < cnst_tri.size(); ++i)
+		for (int c = 0; c < cnst_tri.size(); ++c)
 		{
 			for (t_iter = out_tris.begin(); t_iter != out_tris.end(); )
 			{
 				tmp_tri = *t_iter;
-				if (cnst_tri[i] == tmp_tri)
+				if (cnst_tri[c] == tmp_tri)
 				{
 					t_iter = out_tris.erase(t_iter);
 					break; // no need to search more
@@ -315,10 +321,41 @@ void dem2tin(const std::vector<double> &dem, double xmin, double xmax, double ym
 			}
 		}
 
+		// remove cnst_tri from its circumed DEM's circum triangle list
+		for (int c = 0; c < cnst_tri.size(); ++c)
+		{
+			for (int i = 0; i < cnst_tri[c]->circum_dem.size(); ++i)
+			{
+				tmp_dem = cnst_tri[c]->circum_dem[i];
+				for (t_iter = tmp_dem->circum_host.begin(); t_iter != tmp_dem->circum_host.end(); )
+				{
+					if (cnst_tri[c] == *t_iter)
+					{
+						t_iter = tmp_dem->circum_host.erase(t_iter);
+						break;
+					}
+					else t_iter++;
+				}
+			}
+		}
+
+		// remove dem_grid[0] from its circumed triangle's circum DEM list
+		for (int c = 0; c < cnst_tri.size(); ++c)
+		{
+			for (d_iter = cnst_tri[c]->circum_dem.begin(); d_iter != cnst_tri[c]->circum_dem.end(); )
+			{
+				if (dem_grid[0] == *d_iter)
+				{
+					d_iter = cnst_tri[c]->circum_dem.erase(d_iter);
+					break;
+				}
+				else d_iter++;
+			}
+		}
+
 		// clear host and circumcircle triangles for the used DEM location
-		dem_grid[0]->circum_host.clear();
 		d_iter = dem_grid.begin();
-		tmp_dem = *d_iter; delete tmp_dem;
+		tmp_dem = *d_iter; tmp_dem->circum_host.clear(); delete tmp_dem;
 		dem_grid.erase(d_iter);
 
 		// loop to remove duplicate edges
@@ -360,38 +397,26 @@ void dem2tin(const std::vector<double> &dem, double xmin, double xmax, double ym
 			}
 		}
 
-		// purge circumcircle triangles for all DEM data
+		// loop all DEM data to update host triangles
 		for (int c = 0; c < cnst_tri.size(); ++c)
 		{
-			for (int i = 0; i < dem_grid.size(); ++i)
+			for (int i = 0; i < cnst_tri[c]->circum_dem.size(); ++i)
 			{
-				for (t_iter = dem_grid[i]->circum_host.begin(); t_iter != dem_grid[i]->circum_host.end(); )
+				tmp_dem = cnst_tri[c]->circum_dem[i];
+				for (int n = 0; n < new_tri.size(); ++n) // search in newly created triangles to find new host
 				{
-					if (cnst_tri[c] == *t_iter)
+					if (new_tri[n]->bound_location(tmp_dem->x, tmp_dem->y))
 					{
-						t_iter = dem_grid[i]->circum_host.erase(t_iter);
-						break; // no need to search more
+						tmp_dem->host = new_tri[n];
+						tmp_dem->err = fabs(new_tri[n]->interpolate(tmp_dem->x, tmp_dem->y) - tmp_dem->elev);
+						break; // already found, no need to search more
 					}
-					else t_iter++;
-				}
-			}
-		}
-
-		// loop all DEM data to update host and circumcircle triangles
-		for (int i = 0; i < dem_grid.size(); ++i)
-		{
-			for (int n = 0; n < new_tri.size(); ++n) // search in newly created triangles to find new host
-			{
-				if (new_tri[n]->bound_location(dem_grid[i]->x, dem_grid[i]->y))
-				{
-					dem_grid[i]->host = new_tri[n];
-					dem_grid[i]->err = fabs(new_tri[n]->interpolate(dem_grid[i]->x, dem_grid[i]->y) - dem_grid[i]->elev);
-					break; // already found, no need to search more
 				}
 			}
 		}
 
 		// Find circum_host triangles for all DEM locations
+		// cnst_tri's circum area doesn't over cover new_tri's circum area
 		for (int i = 0; i < dem_grid.size(); ++i)
 		{
 			for (int n = 0; n < new_tri.size(); ++n) // search in newly created triangles to find new circumcircle triangles
@@ -400,16 +425,19 @@ void dem2tin(const std::vector<double> &dem, double xmin, double xmax, double ym
 					+ (new_tri[n]->cy - dem_grid[i]->y) * (new_tri[n]->cy - dem_grid[i]->y);
 				if ((dist - new_tri[n]->cr) <= ZERO) // Points on the circumcircle are also included
 				{
+					new_tri[n]->circum_dem.push_back(dem_grid[i]);
 					dem_grid[i]->circum_host.push_back(new_tri[n]);
 					// no beak here. There might be more than one triangle's circumcircle includes the DEM location
 				}
 			}
 		}
+		
 
 		// destroy memories used by cnst_edge
 		for (int c = 0; c < cnst_tri.size(); ++c)
 		{
 			tmp_tri = cnst_tri[c];
+			tmp_tri->circum_dem.clear();
 			delete tmp_tri; tmp_tri = nullptr;
 		}