diff --git a/stim/biomodels/network.h b/stim/biomodels/network.h
index 06c376e..38eb3a2 100644
--- a/stim/biomodels/network.h
+++ b/stim/biomodels/network.h
@@ -473,6 +473,146 @@ public:
 		return E[e];									//return the specified edge (casting it to a fiber)
 	}
 
+	/// subdivide current network
+	void subdivision() {
+
+		std::vector<unsigned> ori_index;		// original index
+		std::vector<unsigned> new_index;		// new index
+		std::vector<edge> nE;					// new edge
+		std::vector<vertex> nV;					// new vector
+		unsigned id = 0;
+
+		for (unsigned i = 0; i < num_edge; i++) {
+			if (E[i].size() == 2) {				// if current edge can't be subdivided
+				stim::centerline<T> line(2);
+				for (unsigned k = 0; k < 2; k++)
+					line[k] = E[i][k];
+				line.update();
+
+				edge new_edge(line);
+
+				vertex new_vertex = new_edge[0];
+				id = E[i].v[0];
+				auto position = std::find(ori_index.begin(), ori_index.end(), id);
+				if (position == ori_index.end()) {		 // new vertex
+					ori_index.push_back(id);
+					new_index.push_back(nV.size());
+
+					new_vertex.e[0].push_back(nE.size());
+					new_edge.v[0] = nV.size();
+					nV.push_back(new_vertex);			// push back vertex as a new vertex
+				}
+				else {									// existing vertex
+					int k = std::distance(ori_index.begin(), position);
+					new_edge.v[0] = new_index[k];
+					nV[new_index[k]].e[0].push_back(nE.size());
+				}
+
+				new_vertex = new_edge[1];
+				id = E[i].v[1];
+				position = std::find(ori_index.begin(), ori_index.end(), id);
+				if (position == ori_index.end()) {		 // new vertex
+					ori_index.push_back(id);
+					new_index.push_back(nV.size());
+
+					new_vertex.e[1].push_back(nE.size());
+					new_edge.v[1] = nV.size();
+					nV.push_back(new_vertex);			// push back vertex as a new vertex
+				}
+				else {									// existing vertex
+					int k = std::distance(ori_index.begin(), position);
+					new_edge.v[1] = new_index[k];
+					nV[new_index[k]].e[1].push_back(nE.size());
+				}
+
+				nE.push_back(new_edge);
+
+				nE[nE.size() - 1].cylinder<T>::set_r(0, E[i].cylinder<T>::r(0));
+				nE[nE.size() - 1].cylinder<T>::set_r(1, E[i].cylinder<T>::r(1));
+			}
+			else {								// subdivide current edge
+				for (unsigned j = 0; j < E[i].size() - 1; j++) {
+					stim::centerline<T> line(2);
+					for (unsigned k = 0; k < 2; k++)
+						line[k] = E[i][j + k];
+					line.update();
+
+					edge new_edge(line);
+
+					if (j == 0) {						// edge contains original starting point
+						vertex new_vertex = new_edge[0];
+						id = E[i].v[0];
+						auto position = std::find(ori_index.begin(), ori_index.end(), id);
+						if (position == ori_index.end()) {		 // new vertex
+							ori_index.push_back(id);
+							new_index.push_back(nV.size());
+
+							new_vertex.e[0].push_back(nE.size());
+							new_edge.v[0] = nV.size();
+							nV.push_back(new_vertex);			// push back vertex as a new vertex
+						}
+						else {									// existing vertex
+							int k = std::distance(ori_index.begin(), position);
+							new_edge.v[0] = new_index[k];
+							nV[new_index[k]].e[0].push_back(nE.size());
+						}
+
+						new_vertex = new_edge[1];
+						new_vertex.e[1].push_back(nE.size());
+						new_edge.v[1] = nV.size();
+						nV.push_back(new_vertex);				// push back internal point as a new vertex
+
+						nE.push_back(new_edge);
+					}
+
+					else if (j == E[i].size() - 2) {	// edge contains original ending point
+
+						vertex new_vertex = new_edge[1];
+						nV[nV.size() - 1].e[0].push_back(nE.size());
+						new_edge.v[0] = nV.size() - 1;
+
+						id = E[i].v[1];
+						auto position = std::find(ori_index.begin(), ori_index.end(), id);
+						if (position == ori_index.end()) {		 // new vertex
+							ori_index.push_back(id);
+							new_index.push_back(nV.size());
+
+							new_vertex.e[1].push_back(nE.size());
+							new_edge.v[1] = nV.size();
+							nV.push_back(new_vertex);			// push back vertex as a new vertex
+						}
+						else {									// existing vertex
+							int k = std::distance(ori_index.begin(), position);
+							new_edge.v[1] = new_index[k];
+							nV[new_index[k]].e[1].push_back(nE.size());
+						}
+
+						nE.push_back(new_edge);
+					}
+
+					else {
+						vertex new_vertex = new_edge[1];
+
+						nV[nV.size() - 1].e[0].push_back(nE.size());
+						new_vertex.e[1].push_back(nE.size());
+						new_edge.v[0] = nV.size() - 1;
+						new_edge.v[1] = nV.size();
+						nV.push_back(new_vertex);
+
+						nE.push_back(new_edge);
+					}
+
+					// get radii
+					nE[nE.size() - 1].cylinder<T>::set_r(0, E[i].cylinder<T>::r(j));
+					nE[nE.size() - 1].cylinder<T>::set_r(1, E[i].cylinder<T>::r(j + 1));
+				}
+			}
+		}
+
+		(*this).E = nE;
+		(*this).V = nV;
+	}
+
 	//load a network from an OBJ file
 	void load_obj(std::string filename){
 
@@ -715,7 +855,7 @@ public:
 			edge new_edge(C3);								// new edge	
 
 			//create an edge from the given centerline
-			unsigned int I = new_edge.size();				//calculate the number of points on the centerline
+			unsigned int I = (unsigned)new_edge.size();				//calculate the number of points on the centerline
 			
 			//get the first and last vertex IDs for the line
 			i[0] = S.E[l].front();
--
libgit2 0.21.4