diff --git a/FileIO/Legacy/MeshIO.cpp b/FileIO/Legacy/MeshIO.cpp
index aadd737928e275c1d3efaf079ea6a85f392695d4..12e1c4bd2e4f35061df86927b2657b81dc1d36fc 100644
--- a/FileIO/Legacy/MeshIO.cpp
+++ b/FileIO/Legacy/MeshIO.cpp
@@ -173,8 +173,8 @@ MeshLib::Element* MeshIO::readElement(const std::string& line,
 			ss >> idx[i];
 		// edge_nodes array will be deleted from Line object
 		MeshLib::Node** edge_nodes = new MeshLib::Node*[2];
-		edge_nodes[0] = nodes[idx[1]];
-		edge_nodes[1] = nodes[idx[0]];
+		for (unsigned k(0); k < 2; ++k)
+			edge_nodes[k] = nodes[idx[k]];
 		elem = new MeshLib::Line(edge_nodes, patch_index);
 		break;
 	}
@@ -183,7 +183,7 @@ MeshLib::Element* MeshIO::readElement(const std::string& line,
 			ss >> idx[i];
 		MeshLib::Node** tri_nodes = new MeshLib::Node*[3];
 		for (unsigned k(0); k < 3; ++k)
-			tri_nodes[k] = nodes[idx[2 - k]];
+			tri_nodes[k] = nodes[idx[k]];
 		elem = new MeshLib::Tri(tri_nodes, patch_index);
 		break;
 	}
@@ -192,7 +192,7 @@ MeshLib::Element* MeshIO::readElement(const std::string& line,
 			ss >> idx[i];
 		MeshLib::Node** quad_nodes = new MeshLib::Node*[4];
 		for (unsigned k(0); k < 4; ++k)
-			quad_nodes[k] = nodes[idx[3 - k]];
+			quad_nodes[k] = nodes[idx[k]];
 		elem = new MeshLib::Quad(quad_nodes, patch_index);
 		break;
 	}
@@ -201,7 +201,7 @@ MeshLib::Element* MeshIO::readElement(const std::string& line,
 			ss >> idx[i];
 		MeshLib::Node** tet_nodes = new MeshLib::Node*[4];
 		for (unsigned k(0); k < 4; ++k)
-			tet_nodes[k] = nodes[idx[3 - k]];
+			tet_nodes[k] = nodes[idx[k]];
 		elem = new MeshLib::Tet(tet_nodes, patch_index);
 		break;
 	}
@@ -210,7 +210,7 @@ MeshLib::Element* MeshIO::readElement(const std::string& line,
 			ss >> idx[i];
 		MeshLib::Node** hex_nodes = new MeshLib::Node*[8];
 		for (unsigned k(0); k < 8; ++k)
-			hex_nodes[k] = nodes[idx[7 - k]];
+			hex_nodes[k] = nodes[idx[k]];
 		elem = new MeshLib::Hex(hex_nodes, patch_index);
 		break;
 	}
@@ -219,7 +219,7 @@ MeshLib::Element* MeshIO::readElement(const std::string& line,
 			ss >> idx[i];
 		MeshLib::Node** pyramid_nodes = new MeshLib::Node*[5];
 		for (unsigned k(0); k < 5; ++k)
-			pyramid_nodes[k] = nodes[idx[4 - k]];
+			pyramid_nodes[k] = nodes[idx[k]];
 		elem = new MeshLib::Pyramid(pyramid_nodes, patch_index);
 		break;
 	}
@@ -228,12 +228,13 @@ MeshLib::Element* MeshIO::readElement(const std::string& line,
 			ss >> idx[i];
 		MeshLib::Node** prism_nodes = new MeshLib::Node*[6];
 		for (unsigned k(0); k < 6; ++k)
-			prism_nodes[k] = nodes[idx[5 - k]];
+			prism_nodes[k] = nodes[idx[k]];
 		elem = new MeshLib::Prism(prism_nodes, patch_index);
 		break;
 	}
 	default:
 		elem = NULL;
+		break;
 	}
 
 	delete [] idx;
@@ -262,7 +263,7 @@ int MeshIO::write(std::ostream &out)
 	out << "$ELEMENTS\n"
 		<< "  ";
 
-	writeElementsExceptLines(_mesh->getElements(), out);
+	writeElements(_mesh->getElements(), out);
 
 	out << " $LAYER\n"
 		<< "  0\n"
@@ -276,37 +277,18 @@ void MeshIO::setMesh(const MeshLib::Mesh* mesh)
 	_mesh = mesh;
 }
 
-void MeshIO::writeElementsExceptLines(std::vector<MeshLib::Element*> const& ele_vec,
+void MeshIO::writeElements(std::vector<MeshLib::Element*> const& ele_vec,
                                       std::ostream &out)
 {
 	const size_t ele_vector_size (ele_vec.size());
-	const double epsilon (std::numeric_limits<double>::epsilon());
-	std::vector<bool> non_line_element (ele_vector_size, true);
-	std::vector<bool> non_null_element (ele_vector_size, true);
-	size_t n_elements(0);
 
+	out << ele_vector_size << "\n";
 	for (size_t i(0); i < ele_vector_size; ++i) {
-		if ((ele_vec[i])->getGeomType() == MeshElemType::LINE) {
-			non_line_element[i] = false;
-			non_null_element[i] = false;
-		} else {
-			if (ele_vec[i]->getContent() < epsilon) {
-				non_null_element[i] = false;
-			} else {
-				++n_elements;
-			}
-		}
-	}
-	out << n_elements << "\n";
-	for (size_t i(0), k(0); i < ele_vector_size; ++i) {
-		if (non_line_element[i] && non_null_element[i]) {
-			out << k << " " << ele_vec[i]->getValue() << " " << MeshElemType2String(ele_vec[i]->getGeomType()) << " ";
-			unsigned nElemNodes (ele_vec[i]->getNNodes());
-			for(size_t j = 0; j < nElemNodes; ++j)
-				out << ele_vec[i]->getNode(nElemNodes - j - 1)->getID() << " ";
-			out << "\n";
-			++k;
-		}
+		out << i << " " << ele_vec[i]->getValue() << " " << MeshElemType2String(ele_vec[i]->getGeomType()) << " ";
+		unsigned nElemNodes (ele_vec[i]->getNNodes());
+		for(size_t j = 0; j < nElemNodes; ++j)
+			out << ele_vec[i]->getNode(j)->getID() << " ";
+		out << "\n";
 	}
 }
 
diff --git a/FileIO/Legacy/MeshIO.h b/FileIO/Legacy/MeshIO.h
index e8e9d7875a15a36f6f449225d433b22dce87e552..bbcdd6288149736a49866fdc1823276e4183146a 100644
--- a/FileIO/Legacy/MeshIO.h
+++ b/FileIO/Legacy/MeshIO.h
@@ -53,7 +53,7 @@ protected:
 	int write(std::ostream &out);
 
 private:
-	void writeElementsExceptLines(std::vector<MeshLib::Element*> const& ele_vec, std::ostream &out);
+	void writeElements(std::vector<MeshLib::Element*> const& ele_vec, std::ostream &out);
 	MeshLib::Element* readElement(const std::string& line, const std::vector<MeshLib::Node*> &nodes);
 
 	double* _edge_length[2];
diff --git a/Utils/MeshEdit/CMakeLists.txt b/Utils/MeshEdit/CMakeLists.txt
index e637d2cd66809d82561a472f4c9ef31eb3ef9b9a..aae8eb8ce2dee237d637d0dae1ada9f596360482 100644
--- a/Utils/MeshEdit/CMakeLists.txt
+++ b/Utils/MeshEdit/CMakeLists.txt
@@ -33,3 +33,21 @@ IF(QT4_FOUND)
 		PROPERTIES FOLDER Utilities)
 
 ENDIF() # QT4_FOUND
+
+ADD_EXECUTABLE( reverseMeshNodeOrdering reverseMeshNodeOrdering.cpp )
+TARGET_LINK_LIBRARIES( reverseMeshNodeOrdering
+	BaseLib
+	FileIO
+	MathLib
+	MeshLib
+)
+SET_TARGET_PROPERTIES(reverseMeshNodeOrdering PROPERTIES FOLDER Utilities)
+
+ADD_EXECUTABLE( removeMeshElements removeMeshElements.cpp )
+TARGET_LINK_LIBRARIES( removeMeshElements
+	BaseLib
+	FileIO
+	MathLib
+	MeshLib
+)
+SET_TARGET_PROPERTIES(removeMeshElements PROPERTIES FOLDER Utilities)
diff --git a/Utils/MeshEdit/removeMeshElements.cpp b/Utils/MeshEdit/removeMeshElements.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..64abe499358bf68776acc35edb954e8101b377f5
--- /dev/null
+++ b/Utils/MeshEdit/removeMeshElements.cpp
@@ -0,0 +1,192 @@
+/**
+ * @file   removeMeshElements.cpp
+ * @author Norihiro Watanabe
+ * @date   2013/10/15
+ * @brief  Remove mesh elements
+ *
+ * @copyright
+ * Copyright (c) 2013, OpenGeoSys Community (http://www.opengeosys.org)
+ *            Distributed under a Modified BSD License.
+ *              See accompanying file LICENSE.txt or
+ *              http://www.opengeosys.org/LICENSE.txt
+ */
+
+// TCLAP
+#include "tclap/CmdLine.h"
+
+// ThirdParty/logog
+#include "logog/include/logog.hpp"
+
+// BaseLib
+#include "LogogSimpleFormatter.h"
+
+// FileIO
+#include "Legacy/MeshIO.h"
+#include "readMeshFromFile.h"
+
+// MeshLib
+#include "Mesh.h"
+#include "Node.h"
+#include "Elements/Element.h"
+#include "MeshEnums.h"
+
+std::vector<std::size_t> searchByMaterialID(const std::vector<MeshLib::Element*> & ele_vec, unsigned matID)
+{
+	std::vector<std::size_t> matchedIDs;
+	std::size_t i = 0;
+	for (MeshLib::Element* ele : ele_vec) {
+		if (ele->getValue()==matID)
+			matchedIDs.push_back(i);
+		i++;
+	}
+	return matchedIDs;
+}
+
+std::vector<std::size_t> searchByElementType(const std::vector<MeshLib::Element*> & ele_vec, MeshElemType eleType)
+{
+	std::vector<std::size_t> matchedIDs;
+	std::size_t i = 0;
+	for (MeshLib::Element* ele : ele_vec) {
+		if (ele->getGeomType()==eleType)
+			matchedIDs.push_back(i);
+		i++;
+	}
+	return matchedIDs;
+}
+
+std::vector<std::size_t> searchByZeroContent(const std::vector<MeshLib::Element*> & ele_vec)
+{
+	std::vector<std::size_t> matchedIDs;
+	std::size_t i = 0;
+	for (MeshLib::Element* ele : ele_vec) {
+		if (ele->getContent()==.0)
+			matchedIDs.push_back(i);
+		i++;
+	}
+	return matchedIDs;
+}
+
+void updateUnion(const std::vector<std::size_t> &vec1, std::vector<std::size_t> &vec2)
+{
+	std::vector<std::size_t> vec_temp(vec1.size() + vec2.size());
+	auto it = std::set_union(vec1.begin(), vec1.end(), vec2.begin(), vec2.end(), vec_temp.begin());
+	vec_temp.resize(it - vec_temp.begin());
+	vec2.assign(vec_temp.begin(), vec_temp.end());
+}
+
+std::vector<MeshLib::Element*> excludeElements(const std::vector<MeshLib::Element*> & vec_src_eles, const std::vector<std::size_t> &vec_removed)
+{
+	std::vector<MeshLib::Element*> vec_dest_eles(vec_src_eles.size() - vec_removed.size());
+	std::size_t k=0;
+	for (std::size_t i=0; i<vec_src_eles.size(); i++) {
+		if (std::find(vec_removed.begin(), vec_removed.end(), i) == vec_removed.end()) {
+			vec_dest_eles[k] = vec_src_eles[i];
+			k++;
+		}
+	}
+	return vec_dest_eles;
+}
+
+void copyNodesElements(	const std::vector<MeshLib::Node*> src_nodes,
+						const std::vector<MeshLib::Element*> & src_eles,
+						std::vector<MeshLib::Node*> &dst_nodes,
+						std::vector<MeshLib::Element*> &dst_eles)
+{
+	// copy nodes
+	dst_nodes.resize(src_nodes.size());
+	for (std::size_t i=0; i<dst_nodes.size(); i++) {
+		dst_nodes[i] = new MeshLib::Node(*src_nodes[i]);
+	}
+
+	// copy elements with new nodes
+	dst_eles.resize(src_eles.size());
+	for (std::size_t i=0; i<dst_eles.size(); i++) {
+		auto* src_ele = src_eles[i];
+		auto* dst_ele = src_ele->clone();
+		for (unsigned j=0; j<src_ele->getNNodes(); j++) {
+			dst_ele->setNode(j, dst_nodes[src_ele->getNode(j)->getID()]);
+		}
+		dst_eles[i] = dst_ele;
+	}
+}
+
+int main (int argc, char* argv[])
+{
+	LOGOG_INITIALIZE();
+	logog::Cout* logog_cout (new logog::Cout);
+	BaseLib::LogogSimpleFormatter *custom_format (new BaseLib::LogogSimpleFormatter);
+	logog_cout->SetFormatter(*custom_format);
+
+	TCLAP::CmdLine cmd("Remove mesh elements.", ' ', "0.1");
+	TCLAP::ValueArg<std::string> mesh_in("i", "mesh-input-file",
+	                                     "the name of the file containing the input mesh", true,
+	                                     "", "file name of input mesh");
+	cmd.add(mesh_in);
+	TCLAP::ValueArg<std::string> mesh_out("o", "mesh-output-file",
+	                                      "the name of the file the mesh will be written to", true,
+	                                      "", "file name of output mesh");
+	cmd.add(mesh_out);
+	TCLAP::SwitchArg zveArg("z", "zero-volume", "remove zero volume elements", false);
+	cmd.add(zveArg);
+	TCLAP::MultiArg<std::string> eleTypeArg("t", "element-type",
+	                                      "element type to be removed", false, "element type");
+	cmd.add(eleTypeArg);
+	TCLAP::MultiArg<unsigned> matIDArg("m", "material-id",
+	                                      "material id", false, "material id");
+	cmd.add(matIDArg);
+	cmd.parse(argc, argv);
+
+	MeshLib::Mesh* mesh (FileIO::readMeshFromFile(mesh_in.getValue()));
+	INFO("Mesh read: %d nodes, %d elements.", mesh->getNNodes(), mesh->getNElements());
+
+	// search elements IDs to be removed
+	std::vector<std::size_t> vec_elementIDs_removed;
+	if (zveArg.isSet()) {
+		std::vector<std::size_t> vec_matched = searchByZeroContent(mesh->getElements());
+		updateUnion(vec_matched, vec_elementIDs_removed);
+		INFO("%d zero volume elements found.", vec_matched.size());
+	}
+	if (eleTypeArg.isSet()) {
+		std::vector<std::string> eleTypeNames = eleTypeArg.getValue();
+		for (auto typeName : eleTypeNames) {
+			MeshElemType type = String2MeshElemType(typeName);
+			if (type == MeshElemType::INVALID) continue;
+			std::vector<std::size_t> vec_matched = searchByElementType(mesh->getElements(), type);
+			updateUnion(vec_matched, vec_elementIDs_removed);
+			INFO("%d %s elements found.", vec_matched.size(), typeName.c_str());
+		}
+	}
+	if (matIDArg.isSet()) {
+		std::vector<unsigned> vec_matID = matIDArg.getValue();
+		for (auto matID : vec_matID) {
+			std::vector<std::size_t> vec_matched = searchByMaterialID(mesh->getElements(), matID);
+			updateUnion(vec_matched, vec_elementIDs_removed);
+			INFO("%d elements with material ID %d found.", vec_matched.size(), matID);
+		}
+	}
+
+	// remove the elements
+	INFO("Removing total %d elements...", vec_elementIDs_removed.size());
+	std::vector<MeshLib::Element*> tmp_eles = excludeElements(mesh->getElements(), vec_elementIDs_removed);
+	INFO("%d elements remained.", tmp_eles.size());
+	std::vector<MeshLib::Node*> new_nodes;
+	std::vector<MeshLib::Element*> new_eles;
+	copyNodesElements(mesh->getNodes(), tmp_eles, new_nodes, new_eles);
+
+	// create a new mesh object. Unsued nodes are removed while construction
+	MeshLib::Mesh* new_mesh(new MeshLib::Mesh(mesh->getName(), new_nodes, new_eles));
+
+	// write into a file
+	FileIO::MeshIO meshIO;
+	meshIO.setMesh(new_mesh);
+	meshIO.writeToFile(mesh_out.getValue());
+
+	delete custom_format;
+	delete logog_cout;
+	LOGOG_SHUTDOWN();
+
+	return 0;
+}
+
+
+
diff --git a/Utils/MeshEdit/reverseMeshNodeOrdering.cpp b/Utils/MeshEdit/reverseMeshNodeOrdering.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..32b25ba79ac0eb0763de160c7f615fc45102c8fd
--- /dev/null
+++ b/Utils/MeshEdit/reverseMeshNodeOrdering.cpp
@@ -0,0 +1,77 @@
+/**
+ * @file   reverseMeshNodeOrdering.cpp
+ * @author Norihiro Watanabe
+ * @date   2013/10/15
+ * @brief  Reverse element node ordering
+ *
+ * @copyright
+ * Copyright (c) 2013, OpenGeoSys Community (http://www.opengeosys.org)
+ *            Distributed under a Modified BSD License.
+ *              See accompanying file LICENSE.txt or
+ *              http://www.opengeosys.org/LICENSE.txt
+ */
+
+// TCLAP
+#include "tclap/CmdLine.h"
+
+// ThirdParty/logog
+#include "logog/include/logog.hpp"
+
+// BaseLib
+#include "LogogSimpleFormatter.h"
+
+// FileIO
+#include "Legacy/MeshIO.h"
+#include "readMeshFromFile.h"
+
+// MeshLib
+#include "Mesh.h"
+#include "Elements/Element.h"
+
+void reverseNodeOrdering(std::vector<MeshLib::Element*> & ele_vec)
+{
+	for (MeshLib::Element* ele : ele_vec) {
+		unsigned nElemNodes (ele->getNNodes());
+		std::vector<MeshLib::Node*> originalNodes(ele->getNodes(), ele->getNodes() + nElemNodes);
+		for(size_t j = 0; j < nElemNodes; ++j)
+			ele->setNode(j, originalNodes[nElemNodes - j - 1]);
+	}
+}
+
+int main (int argc, char* argv[])
+{
+	LOGOG_INITIALIZE();
+	logog::Cout* logog_cout (new logog::Cout);
+	BaseLib::LogogSimpleFormatter *custom_format (new BaseLib::LogogSimpleFormatter);
+	logog_cout->SetFormatter(*custom_format);
+
+	TCLAP::CmdLine cmd("Reverse the node ordering of mesh elements.", ' ', "0.1");
+	TCLAP::ValueArg<std::string> mesh_in("i", "mesh-input-file",
+	                                     "the name of the file containing the input mesh", true,
+	                                     "", "file name of input mesh");
+	cmd.add(mesh_in);
+	TCLAP::ValueArg<std::string> mesh_out("o", "mesh-output-file",
+	                                      "the name of the file the mesh will be written to", true,
+	                                      "", "file name of output mesh");
+	cmd.add(mesh_out);
+	cmd.parse(argc, argv);
+
+	MeshLib::Mesh* mesh (FileIO::readMeshFromFile(mesh_in.getValue()));
+	INFO("Mesh read: %d nodes, %d elements.", mesh->getNNodes(), mesh->getNElements());
+
+	INFO("Reversing the node ordering...");
+	reverseNodeOrdering(const_cast<std::vector<MeshLib::Element*>&>(mesh->getElements()));
+
+	FileIO::MeshIO meshIO;
+	meshIO.setMesh(mesh);
+	meshIO.writeToFile(mesh_out.getValue());
+
+	delete custom_format;
+	delete logog_cout;
+	LOGOG_SHUTDOWN();
+
+	return 0;
+}
+
+
+