diff --git a/schnapps/plugins/import/CMakeLists.txt b/schnapps/plugins/import/CMakeLists.txt
index 170a89a..79a1381 100644
--- a/schnapps/plugins/import/CMakeLists.txt
+++ b/schnapps/plugins/import/CMakeLists.txt
@@ -8,7 +8,6 @@ find_package(Qt5 COMPONENTS Widgets REQUIRED)
 find_package(tinyxml2 REQUIRED)
 find_package(lm6 REQUIRED)
 find_package(ply REQUIRED)
-find_package(QOGLViewer REQUIRED)
 
 set(HEADER_FILES
 	dll.h
diff --git a/schnapps/plugins/shallow_water/CMakeLists.txt b/schnapps/plugins/shallow_water/CMakeLists.txt
index 34d8e09..5287eb0 100644
--- a/schnapps/plugins/shallow_water/CMakeLists.txt
+++ b/schnapps/plugins/shallow_water/CMakeLists.txt
@@ -37,7 +37,7 @@ target_link_libraries(${PROJECT_NAME}
 	schnapps_core
 	${cgogn_core_LIBRARIES}
 	${cgogn_modeling_LIBRARIES}
-	${cgogn_geoemtry_LIBRARIES}
+	${cgogn_geometry_LIBRARIES}
 	${Qt5Core_LIBRARIES}
 	${Qt5Widgets_LIBRARIES}
 )
diff --git a/schnapps/plugins/shallow_water/shallow_water.cpp b/schnapps/plugins/shallow_water/shallow_water.cpp
index ecf8aa6..16257f3 100644
--- a/schnapps/plugins/shallow_water/shallow_water.cpp
+++ b/schnapps/plugins/shallow_water/shallow_water.cpp
@@ -45,6 +45,7 @@ bool Plugin_ShallowWater::enable()
 
 	map_ = static_cast<CMap2Handler*>(schnapps_->add_map("shallow_water", 2));
 	map2_ = static_cast<CMap2*>(map_->get_map());
+	qtrav_ = cgogn::make_unique<CMap2::QuickTraversor>(*map2_);
 
 	position_ = map_->add_attribute<VEC3, CMap2::Vertex::ORBIT>("position");
 	water_position_ = map_->add_attribute<VEC3, CMap2::Vertex::ORBIT>("water_position");
@@ -55,9 +56,11 @@ bool Plugin_ShallowWater::enable()
 	h_tmp_ = map_->add_attribute<SCALAR, CMap2::Face::ORBIT>("hauteur_tmp");
 	q_ = map_->add_attribute<SCALAR, CMap2::Face::ORBIT>("debit");
 	q_tmp_ = map_->add_attribute<SCALAR, CMap2::Face::ORBIT>("debit_tmp");
+
 	centroid_ = map_->add_attribute<VEC3, CMap2::Face::ORBIT>("centroid");
 	length_ = map_->add_attribute<SCALAR, CMap2::Face::ORBIT>("length");
 	phi_ = map_->add_attribute<SCALAR, CMap2::Face::ORBIT>("phi");
+	zb_ = map_->add_attribute<SCALAR, CMap2::Face::ORBIT>("zb");
 
 	subd_code_ = map_->add_attribute<uint32, CMap2::Face::ORBIT>("subdivision_code");
 
@@ -74,26 +77,52 @@ bool Plugin_ShallowWater::enable()
 	boundaryL_ = CMap2::Edge(grid.vertex_table_[nbc+1].dart);
 	boundaryR_ = CMap2::Edge(grid.vertex_table_[nbc].dart);
 
-	cgogn::geometry::compute_centroid<VEC3, CMap2::Face>(*map2_, position_, centroid_);
+	qtrav_->build<CMap2::Vertex>();
 
-	map2_->parallel_foreach_cell([&] (CMap2::Face f, uint32)
-	{
-		CMap2::Edge e1(f.dart);
-		CMap2::Edge e2(map2_->phi1(f.dart));
-		if (map2_->is_incident_to_boundary(e1))
+	qtrav_->build<CMap2::Edge>([&] (CMap2::Edge e) -> cgogn::Dart {
+		if (map2_->is_incident_to_boundary(e))
+			return map2_->is_boundary(e.dart) ? map2_->phi2(e.dart) : e.dart;
+		else
 		{
-			length_[f] = cgogn::geometry::length<VEC3>(*map2_, e1, position_);
-			phi_[f] = cgogn::geometry::length<VEC3>(*map2_, e2, position_);
+			if (position_[CMap2::Vertex(e.dart)][0] < position_[CMap2::Vertex(map2_->phi_1(e.dart))][0])
+				return e.dart;
+			else
+				return map2_->phi2(e.dart);
+		}
+	});
+
+	qtrav_->build<CMap2::Face>([&] (CMap2::Face f) -> cgogn::Dart {
+		CMap2::Edge e(f.dart);
+		if (map2_->is_incident_to_boundary(e))
+		{
+			if (position_[CMap2::Vertex(f.dart)][0] < position_[CMap2::Vertex(map2_->phi1(f.dart))][0])
+				return map2_->phi_1(f.dart);
+			else
+				return map2_->phi1(f.dart);
 		}
 		else
 		{
-			phi_[f] = cgogn::geometry::length<VEC3>(*map2_, e1, position_);
-			length_[f] = cgogn::geometry::length<VEC3>(*map2_, e2, position_);
+			if (position_[CMap2::Vertex(f.dart)][0] < position_[CMap2::Vertex(map2_->phi_1(f.dart))][0])
+				return f.dart;
+			else
+				return map2_->phi1(map2_->phi1(f.dart));
 		}
-
-		subd_code_[f] = 1;
 	});
 
+	cgogn::geometry::compute_centroid<VEC3, CMap2::Face>(*map2_, position_, centroid_);
+
+	map2_->parallel_foreach_cell(
+		[&] (CMap2::Face f, uint32)
+		{
+			CMap2::Edge e1(f.dart);
+			CMap2::Edge e2(map2_->phi1(f.dart));
+			phi_[f] = cgogn::geometry::length<VEC3>(*map2_, e1, position_);
+			length_[f] = cgogn::geometry::length<VEC3>(*map2_, e2, position_);
+			subd_code_[f] = 1;
+		},
+		*qtrav_
+	);
+
 	init();
 
 	timer_ = new QTimer(this);
@@ -111,41 +140,78 @@ void Plugin_ShallowWater::disable()
 
 void Plugin_ShallowWater::init()
 {
-	map2_->parallel_foreach_cell([&] (CMap2::Face f, uint32)
-	{
-		if (centroid_[f][0] < -75.)
-			h_[f] = 15.;
-		else if (centroid_[f][0] < -50.)
-			h_[f] = 3.;
-		else if (centroid_[f][0] < -25.)
-			h_[f] = 8.;
-		else if (centroid_[f][0] < 0.)
-			h_[f] = 1.;
-		else if (centroid_[f][0] < 25.)
-			h_[f] = 15.;
-		else if (centroid_[f][0] < 50.)
-			h_[f] = 1.;
-		else if (centroid_[f][0] < 75.)
-			h_[f] = 6.;
-		else
-			h_[f] = 1.;
+	map2_->parallel_foreach_cell(
+		[&] (CMap2::Face f, uint32)
+		{
+			if (centroid_[f][0] < -75.)
+			{
+				h_[f] = 10.;
+				zb_[f] = 0.;
+			}
+			else if (centroid_[f][0] < -50.)
+			{
+				h_[f] = 30.;
+				zb_[f] = 0.;
+			}
+			else if (centroid_[f][0] < -25.)
+			{
+				h_[f] = 10.;
+				zb_[f] = 0.;
+			}
+			else if (centroid_[f][0] < 0.)
+			{
+				h_[f] = 10.;
+				zb_[f] = 0.;
+			}
+			else if (centroid_[f][0] < 25.)
+			{
+				h_[f] = 10.;
+				zb_[f] = 0.;
+			}
+			else if (centroid_[f][0] < 50.)
+			{
+				h_[f] = 10.;
+				zb_[f] = 0.;
+			}
+			else if (centroid_[f][0] < 75.)
+			{
+				h_[f] = 10.;
+				zb_[f] = 0.;
+			}
+			else
+			{
+				h_[f] = 10.;
+				zb_[f] = 0.;
+			}
 
-		// initial water flow is 0 for all cells
-		q_[f] = 0.;
-	});
+			// initial water flow is 0 for all cells
+			q_[f] = 0.;
+		},
+		*qtrav_
+	);
 
 	t_ = 0.;
-	dt_ = 0.01;
+	dt_ = 0.1;
 
-	map2_->parallel_foreach_cell([&] (CMap2::Vertex v, uint32)
-	{
-		SCALAR h = 0;
-		uint32 nbf = 0;
-		map2_->foreach_incident_face(v, [&] (CMap2::Face f) { h += h_[f]; ++nbf; });
-		h /= nbf;
-		water_position_[v][2] = h;
-		scalar_value_[v] = h;
-	});
+	map2_->parallel_foreach_cell(
+		[&] (CMap2::Vertex v, uint32)
+		{
+			SCALAR wh = 0, bh = 0;
+			uint32 nbf = 0;
+			map2_->foreach_incident_face(v, [&] (CMap2::Face f)
+			{
+				wh += h_[f];
+				bh += zb_[f];
+				++nbf;
+			});
+			wh /= nbf;
+			bh /= nbf;
+			water_position_[v][2] = wh;
+			position_[v][2] = bh;
+			scalar_value_[v] = wh-bh;
+		},
+		*qtrav_
+	);
 
 	map_->notify_attribute_change(CMap2::Vertex::ORBIT, "scalar_value");
 	map_->notify_attribute_change(CMap2::Vertex::ORBIT, "water_position");
@@ -174,10 +240,31 @@ bool Plugin_ShallowWater::is_running()
 	return timer_->isActive();
 }
 
+void Plugin_ShallowWater::update_time_step()
+{
+	std::vector<SCALAR> min_dt_per_thread(cgogn::nb_threads());
+	for (SCALAR& d : min_dt_per_thread) d = 0.1;
+
+	map2_->parallel_foreach_cell(
+		[&] (CMap2::Face f, uint32 idx)
+		{
+			SCALAR swept = std::fabs(q_[f]) / std::max(h_[f], 1e-10) + std::sqrt(9.81 * h_[f]);
+			double dt = length_[f] / std::max(swept, 1e-10);
+			min_dt_per_thread[idx] = dt < min_dt_per_thread[idx] ? dt : min_dt_per_thread[idx];
+		},
+		*qtrav_
+	);
+
+	dt_ = *(std::min_element(min_dt_per_thread.begin(), min_dt_per_thread.end()));
+}
+
 void Plugin_ShallowWater::execute_time_step()
 {
 	connectivity_changed_ = false;
 
+	update_time_step();
+//	std::cout << dt_ << std::endl;
+
 	f1_[boundaryL_] = 0.;
 	f2_[boundaryL_] = 5e-1 * 9.81 * phi_[CMap2::Face(boundaryL_.dart)] * (h_[CMap2::Face(boundaryL_.dart)] * h_[CMap2::Face(boundaryL_.dart)]);
 	s0L_[boundaryL_] = 0.;
@@ -186,16 +273,19 @@ void Plugin_ShallowWater::execute_time_step()
 	map2_->parallel_foreach_cell(
 		[&] (CMap2::Edge e, uint32)
 		{
-			std::pair<CMap2::Face, CMap2::Face> faces = get_LR_faces(e);
-			CMap2::Face fL = faces.first;
-			CMap2::Face fR = faces.second;
-			struct Flux F = Solv_HLL(0., 0., phi_[fL], phi_[fR], h_[fL], h_[fR], q_[fL], q_[fR], 1e-3, 9.81);
+			if (map2_->is_incident_to_boundary(e))
+				return;
+
+			CMap2::Face fL(map2_->phi2(e.dart));
+			CMap2::Face fR(e.dart);
+
+			struct Flux F = Solv_HLL(zb_[fL], zb_[fR], phi_[fL], phi_[fR], h_[fL], h_[fR], q_[fL], q_[fR], 1e-3, 9.81);
 			f1_[e] = F.F1;
 			f2_[e] = F.F2;
 			s0L_[e] = F.S0L;
 			s0R_[e] = F.S0R;
 		},
-		[&] (CMap2::Edge e) { return !map2_->is_incident_to_boundary(e); }
+		*qtrav_
 	);
 
 	f1_[boundaryR_] = 0.;
@@ -203,14 +293,16 @@ void Plugin_ShallowWater::execute_time_step()
 	s0L_[boundaryR_] = 0.;
 	s0R_[boundaryR_] = 0.;
 
-	map2_->parallel_foreach_cell([&] (CMap2::Face f, uint32)
-	{
-		std::pair<CMap2::Edge, CMap2::Edge> edges = get_LR_edges(f);
-		CMap2::Edge eL = edges.first;
-		CMap2::Edge eR = edges.second;
-		h_tmp_[f] = h_[f] + (dt_ / (length_[f] * phi_[f])) * (f1_[eL] - f1_[eR]);
-		q_tmp_[f] = q_[f] + (dt_ / (length_[f] * phi_[f])) * ((f2_[eL] - s0L_[eL]) - (f2_[eR] - s0R_[eR]));
-	});
+	map2_->parallel_foreach_cell(
+		[&] (CMap2::Face f, uint32)
+		{
+			CMap2::Edge eL(f.dart);
+			CMap2::Edge eR(map2_->phi1(map2_->phi1(f.dart)));
+			h_tmp_[f] = h_[f] + (dt_ / (length_[f] * phi_[f])) * (f1_[eL] - f1_[eR]);
+			q_tmp_[f] = q_[f] + (dt_ / (length_[f] * phi_[f])) * ((f2_[eL] - s0L_[eL]) - (f2_[eR] - s0R_[eR]));
+		},
+		*qtrav_
+	);
 
 	map2_->swap_attributes(h_, h_tmp_);
 	map2_->swap_attributes(q_, q_tmp_);
@@ -218,15 +310,25 @@ void Plugin_ShallowWater::execute_time_step()
 	try_simplification();
 	try_subdivision();
 
-	map2_->parallel_foreach_cell([&] (CMap2::Vertex v, uint32)
-	{
-		SCALAR h = 0;
-		uint32 nbf = 0;
-		map2_->foreach_incident_face(v, [&] (CMap2::Face f) { h += h_[f]; ++nbf; });
-		h /= nbf;
-		water_position_[v][2] = h;
-		scalar_value_[v] = h;
-	});
+	map2_->parallel_foreach_cell(
+		[&] (CMap2::Vertex v, uint32)
+		{
+			SCALAR wh = 0, bh = 0;
+			uint32 nbf = 0;
+			map2_->foreach_incident_face(v, [&] (CMap2::Face f)
+			{
+				wh += h_[f];
+				bh += zb_[f];
+				++nbf;
+			});
+			wh /= nbf;
+			bh /= nbf;
+			water_position_[v][2] = wh;
+			position_[v][2] = bh;
+			scalar_value_[v] = wh-bh;
+		},
+		*qtrav_
+	);
 
 	if (connectivity_changed_)
 		map_->notify_connectivity_change();
@@ -240,19 +342,18 @@ void Plugin_ShallowWater::execute_time_step()
 void Plugin_ShallowWater::try_subdivision()
 {
 	CMap2::CellCache cache(*map2_);
-	cache.build<CMap2::Face>();
+	cache.build<CMap2::Face>(*qtrav_);
 
 	map2_->foreach_cell(
 		[&] (CMap2::Face f)
 		{
-			std::pair<CMap2::Edge, CMap2::Edge> edges = get_LR_edges(f);
-			CMap2::Edge eL = edges.first;
-			CMap2::Edge eR = edges.second;
+			CMap2::Edge eL(f.dart);
+			CMap2::Edge eR(map2_->phi1(map2_->phi1(f.dart)));
 
 			if (!map2_->is_incident_to_boundary(eL) && !map2_->is_incident_to_boundary(eR))
 			{
 				SCALAR g = std::abs(h_[CMap2::Face(map2_->phi2(eL.dart))] - h_[CMap2::Face(map2_->phi2(eR.dart))]);
-				if (g > 1. && subd_code_[f] < (1 << 3))
+				if (g > 1. && subd_code_[f] < (1 << 4))
 				{
 					subdivide_face(f);
 					connectivity_changed_ = true;
@@ -268,9 +369,11 @@ void Plugin_ShallowWater::try_simplification()
 	map2_->foreach_cell(
 		[&] (CMap2::Face f)
 		{
-			std::pair<CMap2::Edge, CMap2::Edge> edges = get_LR_edges(f);
-			CMap2::Edge eL = edges.first;
-			CMap2::Edge eR = edges.second;
+			if (subd_code_[f] % 2 != 0)
+				return;
+
+			CMap2::Edge eL(f.dart);
+			CMap2::Edge eR(map2_->phi1(map2_->phi1(f.dart)));
 
 			if (!map2_->is_incident_to_boundary(eR))
 			{
@@ -286,7 +389,7 @@ void Plugin_ShallowWater::try_simplification()
 				}
 			}
 		},
-		[&] (CMap2::Face f) { return subd_code_[f] % 2 == 0; }
+		*qtrav_
 	);
 }
 
@@ -296,17 +399,26 @@ void Plugin_ShallowWater::subdivide_face(CMap2::Face f)
 
 	SCALAR old_h = h_[f];
 	SCALAR old_q = q_[f];
+	SCALAR old_zb = zb_[f];
 
-	std::pair<CMap2::Edge, CMap2::Edge> edges = get_LR_edges(f);
-	CMap2::Edge eL = edges.first;
-	CMap2::Edge eR = edges.second;
+	CMap2::Edge eL(f.dart);
+	CMap2::Edge eR(map2_->phi1(map2_->phi1(f.dart)));
 
 	CMap2::Vertex v1 = map2_->cut_edge(CMap2::Edge(map2_->phi1(eL.dart)));
 	CMap2::Vertex v2 = map2_->cut_edge(CMap2::Edge(map2_->phi1(eR.dart)));
-	CMap2::Edge e = map2_->cut_face(v1.dart, v2.dart);
+	CMap2::Edge e = map2_->cut_face(v2.dart, v1.dart);
+
+	qtrav_->update(v1);
+	qtrav_->update(v2);
+
+	qtrav_->update(CMap2::Edge(map2_->phi2(e.dart)));
+	qtrav_->update(CMap2::Edge(map2_->phi1(e.dart)));
+	qtrav_->update(CMap2::Edge(map2_->phi1(map2_->phi2(e.dart))));
 
 	CMap2::Face fL(eL.dart);
-	CMap2::Face fR(eR.dart);
+	CMap2::Face fR(map2_->phi2(e.dart));
+
+	qtrav_->update(fR);
 
 	subd_code_[fL] = old_subd_code * 2;
 	subd_code_[fR] = old_subd_code * 2 + 1;
@@ -331,6 +443,9 @@ void Plugin_ShallowWater::subdivide_face(CMap2::Face f)
 
 	phi_[fL] = cgogn::geometry::length<VEC3>(*map2_, eL, position_);
 	phi_[fR] = cgogn::geometry::length<VEC3>(*map2_, eR, position_);
+
+	zb_[fL] = old_zb;
+	zb_[fR] = old_zb;
 }
 
 void Plugin_ShallowWater::remove_edge(CMap2::Edge e)
@@ -345,6 +460,8 @@ void Plugin_ShallowWater::remove_edge(CMap2::Edge e)
 	SCALAR old_h_2 = h_[f2];
 	SCALAR old_q_1 = q_[f1];
 	SCALAR old_q_2 = q_[f2];
+	SCALAR old_zb_1 = zb_[f1];
+	SCALAR old_zb_2 = zb_[f2];
 
 	cgogn::Dart d1 = map2_->phi_1(e.dart);
 	cgogn::Dart d2 = map2_->phi_1(map2_->phi2(e.dart));
@@ -365,75 +482,76 @@ void Plugin_ShallowWater::remove_edge(CMap2::Edge e)
 
 	h_[f] = old_h_1;
 	q_[f] = old_q_1;
+	zb_[f] = old_zb_1;
 
 	centroid_[f] = cgogn::geometry::centroid<VEC3>(*map2_, f, position_);
 	length_[f] = cgogn::geometry::length<VEC3>(*map2_, CMap2::Edge(d1), position_);
 	phi_[f] = cgogn::geometry::length<VEC3>(*map2_, CMap2::Edge(map2_->phi_1(d1)), position_);
 }
 
-std::pair<CMap2::Edge, CMap2::Edge> Plugin_ShallowWater::get_LR_edges(CMap2::Face f)
-{
-	CMap2::Edge eL, eR;
-	CMap2::Edge e(f.dart);
-	if (map2_->is_incident_to_boundary(e))
-	{
-		if (map2_->same_cell(e, boundaryL_))
-		{
-			eL = e;
-			eR = CMap2::Edge(map2_->phi1(map2_->phi1(eL.dart)));
-		}
-		else if (map2_->same_cell(e, boundaryR_))
-		{
-			eR = e;
-			eL = CMap2::Edge(map2_->phi1(map2_->phi1(eR.dart)));
-		}
-		else
-		{
-			if (position_[CMap2::Vertex(f.dart)][0] < position_[CMap2::Vertex(map2_->phi1(f.dart))][0])
-			{
-				eL = CMap2::Edge(map2_->phi_1(f.dart));
-				eR = CMap2::Edge(map2_->phi1(f.dart));
-			}
-			else
-			{
-				eL = CMap2::Edge(map2_->phi1(f.dart));
-				eR = CMap2::Edge(map2_->phi_1(f.dart));
-			}
-		}
-	}
-	else
-	{
-		if (position_[CMap2::Vertex(f.dart)][0] < position_[CMap2::Vertex(map2_->phi_1(f.dart))][0])
-		{
-			eL = CMap2::Edge(f.dart);
-			eR = CMap2::Edge(map2_->phi1(map2_->phi1(f.dart)));
-		}
-		else
-		{
-			eL = CMap2::Edge(map2_->phi1(map2_->phi1(f.dart)));
-			eR = CMap2::Edge(f.dart);
-		}
-	}
-
-	return std::make_pair(eL, eR);
-}
-
-std::pair<CMap2::Face, CMap2::Face> Plugin_ShallowWater::get_LR_faces(CMap2::Edge e)
-{
-	CMap2::Face fL, fR;
-	if (position_[CMap2::Vertex(e.dart)][0] < position_[CMap2::Vertex(map2_->phi_1(e.dart))][0])
-	{
-		fL = CMap2::Face(map2_->phi2(e.dart));
-		fR = CMap2::Face(e.dart);
-	}
-	else
-	{
-		fL = CMap2::Face(e.dart);
-		fR = CMap2::Face(map2_->phi2(e.dart));
-	}
-
-	return std::make_pair(fL, fR);
-}
+//std::pair<CMap2::Edge, CMap2::Edge> Plugin_ShallowWater::get_LR_edges(CMap2::Face f)
+//{
+//	CMap2::Edge eL, eR;
+//	CMap2::Edge e(f.dart);
+//	if (map2_->is_incident_to_boundary(e))
+//	{
+//		if (map2_->same_cell(e, boundaryL_))
+//		{
+//			eL = e;
+//			eR = CMap2::Edge(map2_->phi1(map2_->phi1(eL.dart)));
+//		}
+//		else if (map2_->same_cell(e, boundaryR_))
+//		{
+//			eR = e;
+//			eL = CMap2::Edge(map2_->phi1(map2_->phi1(eR.dart)));
+//		}
+//		else
+//		{
+//			if (position_[CMap2::Vertex(f.dart)][0] < position_[CMap2::Vertex(map2_->phi1(f.dart))][0])
+//			{
+//				eL = CMap2::Edge(map2_->phi_1(f.dart));
+//				eR = CMap2::Edge(map2_->phi1(f.dart));
+//			}
+//			else
+//			{
+//				eL = CMap2::Edge(map2_->phi1(f.dart));
+//				eR = CMap2::Edge(map2_->phi_1(f.dart));
+//			}
+//		}
+//	}
+//	else
+//	{
+//		if (position_[CMap2::Vertex(f.dart)][0] < position_[CMap2::Vertex(map2_->phi_1(f.dart))][0])
+//		{
+//			eL = CMap2::Edge(f.dart);
+//			eR = CMap2::Edge(map2_->phi1(map2_->phi1(f.dart)));
+//		}
+//		else
+//		{
+//			eL = CMap2::Edge(map2_->phi1(map2_->phi1(f.dart)));
+//			eR = CMap2::Edge(f.dart);
+//		}
+//	}
+
+//	return std::make_pair(eL, eR);
+//}
+
+//std::pair<CMap2::Face, CMap2::Face> Plugin_ShallowWater::get_LR_faces(CMap2::Edge e)
+//{
+//	CMap2::Face fL, fR;
+//	if (position_[CMap2::Vertex(e.dart)][0] < position_[CMap2::Vertex(map2_->phi_1(e.dart))][0])
+//	{
+//		fL = CMap2::Face(map2_->phi2(e.dart));
+//		fR = CMap2::Face(e.dart);
+//	}
+//	else
+//	{
+//		fL = CMap2::Face(e.dart);
+//		fR = CMap2::Face(map2_->phi2(e.dart));
+//	}
+
+//	return std::make_pair(fL, fR);
+//}
 
 struct Plugin_ShallowWater::Flux Plugin_ShallowWater::Solv_HLL(
 	SCALAR zbL, SCALAR zbR,
diff --git a/schnapps/plugins/shallow_water/shallow_water.h b/schnapps/plugins/shallow_water/shallow_water.h
index 6c5f79b..fe566fc 100644
--- a/schnapps/plugins/shallow_water/shallow_water.h
+++ b/schnapps/plugins/shallow_water/shallow_water.h
@@ -64,6 +64,7 @@ public slots:
 
 private slots:
 
+	void update_time_step();
 	void execute_time_step();
 	void try_subdivision();
 	void try_simplification();
@@ -88,8 +89,8 @@ private slots:
 		SCALAR hmin, SCALAR g
 	);
 
-	std::pair<CMap2::Edge, CMap2::Edge> get_LR_edges(CMap2::Face f);
-	std::pair<CMap2::Face, CMap2::Face> get_LR_faces(CMap2::Edge e);
+//	std::pair<CMap2::Edge, CMap2::Edge> get_LR_edges(CMap2::Face f);
+//	std::pair<CMap2::Face, CMap2::Face> get_LR_faces(CMap2::Edge e);
 
 	ShallowWater_DockTab* dock_tab_;
 
@@ -101,18 +102,22 @@ private slots:
 	CMap2Handler* map_;
 	CMap2* map2_;
 	CMap2::Edge boundaryL_, boundaryR_;
+	std::unique_ptr<CMap2::QuickTraversor> qtrav_;
 
 	CMap2::VertexAttribute<VEC3> position_; // vertices position
 	CMap2::VertexAttribute<VEC3> water_position_;
+
 	CMap2::VertexAttribute<SCALAR> scalar_value_;
 
 	CMap2::FaceAttribute<SCALAR> h_;        // water height
 	CMap2::FaceAttribute<SCALAR> h_tmp_;
 	CMap2::FaceAttribute<SCALAR> q_;        // water flow
 	CMap2::FaceAttribute<SCALAR> q_tmp_;
+
 	CMap2::FaceAttribute<VEC3> centroid_;   // cell centroid
 	CMap2::FaceAttribute<SCALAR> length_;   // cell length
 	CMap2::FaceAttribute<SCALAR> phi_;      // cell width
+	CMap2::FaceAttribute<SCALAR> zb_;       // cell bottom height
 
 	CMap2::FaceAttribute<uint32> subd_code_;// subdivision code
 
diff --git a/schnapps/plugins/shallow_water_2/CMakeLists.txt b/schnapps/plugins/shallow_water_2/CMakeLists.txt
index 82765f7..dfc7411 100644
--- a/schnapps/plugins/shallow_water_2/CMakeLists.txt
+++ b/schnapps/plugins/shallow_water_2/CMakeLists.txt
@@ -3,6 +3,7 @@ project(plugin_shallow_water_2
 )
 
 find_package(cgogn_core REQUIRED)
+find_package(cgogn_io REQUIRED)
 find_package(cgogn_modeling REQUIRED)
 find_package(cgogn_geometry REQUIRED)
 find_package(Qt5 COMPONENTS Core Widgets REQUIRED)
@@ -36,8 +37,9 @@ target_include_directories(${PROJECT_NAME} PUBLIC
 target_link_libraries(${PROJECT_NAME}
 	schnapps_core
 	${cgogn_core_LIBRARIES}
+	${cgogn_io_LIBRARIES}
 	${cgogn_modeling_LIBRARIES}
-	${cgogn_geoemtry_LIBRARIES}
+	${cgogn_geometry_LIBRARIES}
 	${Qt5Core_LIBRARIES}
 	${Qt5Widgets_LIBRARIES}
 )
diff --git a/schnapps/plugins/shallow_water_2/shallow_water.cpp b/schnapps/plugins/shallow_water_2/shallow_water.cpp
index cc258b3..3836304 100644
--- a/schnapps/plugins/shallow_water_2/shallow_water.cpp
+++ b/schnapps/plugins/shallow_water_2/shallow_water.cpp
@@ -27,9 +27,12 @@
 #include <schnapps/core/camera.h>
 
 #include <cgogn/modeling/tiling/triangular_grid.h>
+
 #include <cgogn/geometry/algos/centroid.h>
 #include <cgogn/geometry/algos/length.h>
 
+#include <cgogn/io/map_import.h>
+
 namespace schnapps
 {
 
@@ -46,19 +49,41 @@ bool Plugin_ShallowWater::enable()
 	map_ = static_cast<CMap2Handler*>(schnapps_->add_map("shallow_water_2", 2));
 	map2_ = static_cast<CMap2*>(map_->get_map());
 
-	position_ = map_->add_attribute<VEC3, CMap2::Vertex::ORBIT>("position");
+	cgogn::io::import_surface<VEC3>(*map2_, "/home/kraemer/Desktop/2D.2dm");
+
+	position_ = map_->get_attribute<VEC3, CMap2::Vertex::ORBIT>("position");
+	for (VEC3& p : position_)
+	{
+		p[2] *= 100.;
+	}
 
 	dart_level_ = map_->add_attribute<uint8, CMap2::CDart::ORBIT>("dart_level");
 	face_subd_id_ = map_->add_attribute<uint32, CMap2::Face::ORBIT>("face_subdivision_id");
+	tri_face_ = map_->add_attribute<bool, CMap2::Face::ORBIT>("tri_face");
 
-	cgogn::modeling::TriangularGrid<CMap2> grid(*map2_, nbc, nbc);
-	grid.embed_into_grid(position_, 200.0f, 200.0f, 0.0f);
+//	cgogn::modeling::TriangularGrid<CMap2> grid(*map2_, nbc, nbc);
+//	grid.embed_into_grid(position_, 200.0f, 200.0f, 0.0f);
 
-	map2_->parallel_foreach_cell([&] (CMap2::Face f, uint32)
+	CMap2::CellCache cache(*map2_);
+	cache.build<CMap2::Face>();
+	map2_->foreach_cell([&] (CMap2::Face f)
 	{
+		if (map2_->codegree(f) == 3)
+		{
+			tri_face_[f] = true;
+		}
+		else
+		{
+			tri_face_[f] = false;
+		}
+//		cgogn::Dart d = map2_->phi<11>(f.dart);
+//		map2_->cut_face(f.dart, d);
 		face_subd_id_[f] = 0;
+//		face_subd_id_[CMap2::Face(d)] = 0;
 		map2_->foreach_dart_of_orbit(f, [&] (cgogn::Dart d) { dart_level_[d] = 0; });
-	});
+//		map2_->foreach_dart_of_orbit(CMap2::Face(d), [&] (cgogn::Dart d) { dart_level_[d] = 0; });
+	},
+	cache);
 
 	init();
 
@@ -164,17 +189,16 @@ void Plugin_ShallowWater::try_simplification()
 				connectivity_changed_ = true;
 			}
 		},
-		[&] (CMap2::Face f) -> bool // check connectivity condition
+		[&] (CMap2::Face f) -> bool // check simplification condition
 		{
 			if (face_level(f) == 0)
 				return false;
 
-			if (map2_->codegree(f) != 3)
-				return false;
-
 			switch (face_type(f))
 			{
-				case CORNER: {
+				case TRI_CORNER: {
+					if (map2_->codegree(f) != 3)
+						return false;
 					cgogn::Dart it = oldest_dart(f);
 					it = map2_->phi<12>(it); // central face
 					if (map2_->codegree(CMap2::Face(it)) != 3)
@@ -187,7 +211,9 @@ void Plugin_ShallowWater::try_simplification()
 						return false;
 					break;
 				}
-				case CENTRAL: {
+				case TRI_CENTRAL: {
+					if (map2_->codegree(f) != 3)
+						return false;
 					cgogn::Dart it = f.dart;
 					if (map2_->codegree(CMap2::Face(map2_->phi2(it))) != 3) // corner face 1
 						return false;
@@ -199,6 +225,10 @@ void Plugin_ShallowWater::try_simplification()
 						return false;
 					break;
 				}
+				case QUAD: {
+					return false;
+					break;
+				}
 			}
 
 			return true;
@@ -235,26 +265,37 @@ void Plugin_ShallowWater::subdivide_face(CMap2::Face f)
 		it = next;
 	} while (it != f.dart);
 
-	// cut face into 4 triangles
-	it = map2_->phi1(it);
-	cgogn::Dart it2 = map2_->template phi<11>(it);
-	CMap2::Edge e = map2_->cut_face(it, it2);
-	dart_level_[e.dart] = fl+1;
-	dart_level_[map2_->phi2(e.dart)] = fl+1;
-	face_subd_id_[CMap2::Face(it)] = 4*fid+1;
-
-	it = map2_->template phi<11>(it2);
-	e = map2_->cut_face(it, it2);
-	dart_level_[e.dart] = fl+1;
-	dart_level_[map2_->phi2(e.dart)] = fl+1;
-	face_subd_id_[CMap2::Face(it2)] = 4*fid+2;
-
-	it2 = map2_->template phi<11>(it);
-	e = map2_->cut_face(it, it2);
-	dart_level_[e.dart] = fl+1;
-	dart_level_[map2_->phi2(e.dart)] = fl+1;
-	face_subd_id_[CMap2::Face(it)] = 4*fid+3;
-	face_subd_id_[CMap2::Face(it2)] = 4*fid+4;
+	if (tri_face_[f])
+	{
+		// cut face into 4 triangles
+		it = map2_->phi1(it);
+		cgogn::Dart it2 = map2_->template phi<11>(it);
+		CMap2::Edge e = map2_->cut_face(it, it2);
+		dart_level_[e.dart] = fl+1;
+		dart_level_[map2_->phi2(e.dart)] = fl+1;
+		face_subd_id_[CMap2::Face(it)] = 4*fid+1;
+		tri_face_[CMap2::Face(it)] = true;
+
+		it = map2_->template phi<11>(it2);
+		e = map2_->cut_face(it, it2);
+		dart_level_[e.dart] = fl+1;
+		dart_level_[map2_->phi2(e.dart)] = fl+1;
+		face_subd_id_[CMap2::Face(it2)] = 4*fid+2;
+		tri_face_[CMap2::Face(it2)] = true;
+
+		it2 = map2_->template phi<11>(it);
+		e = map2_->cut_face(it, it2);
+		dart_level_[e.dart] = fl+1;
+		dart_level_[map2_->phi2(e.dart)] = fl+1;
+		face_subd_id_[CMap2::Face(it)] = 4*fid+3;
+		tri_face_[CMap2::Face(it)] = true;
+		face_subd_id_[CMap2::Face(it2)] = 4*fid+4;
+		tri_face_[CMap2::Face(it2)] = true;
+	}
+	else
+	{
+		// cut face into 4 quads
+	}
 }
 
 void Plugin_ShallowWater::simplify_face(CMap2::Face f)
@@ -268,13 +309,13 @@ void Plugin_ShallowWater::simplify_face(CMap2::Face f)
 	cgogn::Dart it = f.dart;
 	switch (face_type(f))
 	{
-		case CORNER: {
+		case TRI_CORNER: {
 			cgogn::Dart od = oldest_dart(f);
 			it = map2_->phi<12>(od); // put 'it' in the central face
 			resF = od;
 			break;
 		}
-		case CENTRAL:
+		case TRI_CENTRAL:
 			resF = map2_->phi_1(map2_->phi2(f.dart));
 			break;
 	}
@@ -326,10 +367,12 @@ uint8 Plugin_ShallowWater::face_level(CMap2::Face f)
 
 Plugin_ShallowWater::FaceType Plugin_ShallowWater::face_type(CMap2::Face f)
 {
-	if (face_subd_id_[f] % 4 == 0)
-		return CENTRAL;
+	if (!tri_face_[f])
+		return QUAD;
+	else if (face_subd_id_[f] % 4 == 0)
+		return TRI_CENTRAL;
 	else
-		return CORNER;
+		return TRI_CORNER;
 }
 
 } // namespace plugin_shallow_water_2
diff --git a/schnapps/plugins/shallow_water_2/shallow_water.h b/schnapps/plugins/shallow_water_2/shallow_water.h
index 26e20aa..56275fb 100644
--- a/schnapps/plugins/shallow_water_2/shallow_water.h
+++ b/schnapps/plugins/shallow_water_2/shallow_water.h
@@ -70,8 +70,9 @@ private slots:
 
 	enum FaceType: uint8
 	{
-		CORNER = 0,
-		CENTRAL
+		TRI_CORNER = 0,
+		TRI_CENTRAL,
+		QUAD
 	};
 
 	void try_subdivision();
@@ -98,6 +99,7 @@ private slots:
 
 	CMap2::CDartAttribute<uint8> dart_level_; // dart insertion level
 	CMap2::FaceAttribute<uint32> face_subd_id_; // face subdivision id
+	CMap2::FaceAttribute<bool> tri_face_; // face is triangle or not
 };
 
 } // namespace plugin_shallow_water_2