diff --git a/.travis.yml b/.travis.yml
index 3e537069b5d7e685009835e0c77b963729d3e2d1..292d2cad5881e6e3d1031c3a4ebb958c71e0b2a0 100644
--- a/.travis.yml
+++ b/.travis.yml
@@ -1,30 +1,46 @@
 language: cpp
+
 compiler:
   - gcc
   - clang
+
 branches:
   only:
     - master
+
 env:
   - CASE=CLI CMAKE_ARGS="-DOGS_BUILD_GUI=OFF -DOGS_BUILD_UTILS=ON"
   - CASE=CLI_PETSC CMAKE_ARGS="-DOGS_BUILD_GUI=OFF -DOGS_BUILD_UTILS=OFF -DOGS_USE_PETSC=ON -DPETSC_DIR=/usr/lib/petscdir/3.4.2/"
   - CASE=GUI CMAKE_ARGS="-DOGS_BUILD_GUI=ON -DVTK_DIR=/home/travis/build/ufz/ogs/VTK-Install/lib/cmake/vtk-6.1/"
+
 before_install:
-  - travis_retry sudo apt-get update; travis_retry sudo apt-get install libeigen3-dev
-  - if [[ "$CC" =~ "gcc" ]]; then travis_retry sudo apt-get install libboost-dev libboost-date-time-dev libboost-filesystem-dev libboost-program-options-dev libboost-system-dev; fi
-  - "wget https://launchpad.net/ubuntu/+source/cmake/2.8.8-2ubuntu1/+build/3441442/+files/cmake_2.8.8-2ubuntu1_amd64.deb"
-  - "wget https://launchpad.net/ubuntu/+archive/primary/+files/cmake-data_2.8.8-2ubuntu1_all.deb"
-  - "sudo apt-get remove cmake-data cmake"
-  - "sudo dpkg --install cmake-data_2.8.8-2ubuntu1_all.deb cmake_2.8.8-2ubuntu1_amd64.deb"
+  # -- External package sources --
+  - sudo add-apt-repository --yes ppa:boost-latest
+  - if [[ "$CASE" == "CLI_PETSC" ]]; then sudo add-apt-repository --yes ppa:fenics-packages/fenics-dev; fi
+  - travis_retry sudo apt-get update;
+
+  # -- Install packages --
+  - travis_retry sudo apt-get install libeigen3-dev
+
+  # Boost
+  - travis_retry sudo apt-get install libboost1.55-dev libboost-date-time1.55-dev libboost-filesystem1.55-dev libboost-program-options1.55-dev libboost-system1.55-dev
+
+  # CMake
+  - travis_retry wget https://launchpad.net/ubuntu/+source/cmake/2.8.8-2ubuntu1/+build/3441442/+files/cmake_2.8.8-2ubuntu1_amd64.deb; travis_retry wget https://launchpad.net/ubuntu/+archive/primary/+files/cmake-data_2.8.8-2ubuntu1_all.deb
+  - sudo apt-get remove cmake-data cmake; sudo dpkg --install cmake-data_2.8.8-2ubuntu1_all.deb cmake_2.8.8-2ubuntu1_amd64.deb
+
+  # Qt and VTK
   - if [[ "$CASE" == "GUI" ]]; then travis_retry sudo apt-get install qt4-dev-tools libshp-dev libgeotiff-dev libxt-dev; fi
-  - if [[ "$CASE" == "GUI" ]]; then wget http://www.opengeosys.org/images/dev/vtk-6.1.0.tar.gz; fi
+  - if [[ "$CASE" == "GUI" ]]; then travis_retry wget http://www.opengeosys.org/images/dev/vtk-6.1.0.tar.gz; fi
   - if [[ "$CASE" == "GUI" ]]; then tar -xf vtk-6.1.0.tar.gz; fi
-  - if [[ "$CASE" == "CLI_PETSC" ]]; then sudo add-apt-repository --yes ppa:amcg/petsc3.4; fi
-  - if [[ "$CASE" == "CLI_PETSC" ]]; then sudo apt-get update; fi
-  - if [[ "$CASE" == "CLI_PETSC" ]]; then sudo apt-get install libhdf5-openmpi-7 libpetsc3.4.2 libpetsc3.4.2-dev; fi
+
+  # PetSc
+  - if [[ "$CASE" == "CLI_PETSC" ]]; then travis_retry sudo apt-get install libpetsc3.4.2 libpetsc3.4.2-dev; fi
+
 script:
   - "pwd & mkdir build && cd build && cmake $CMAKE_ARGS .. && cmake .. && make"
   - make test
+
 notifications:
   email:
     - lars.bilke@ufz.de
diff --git a/FileIO/XmlIO/Qt/XMLQtInterface.cpp b/FileIO/XmlIO/Qt/XMLQtInterface.cpp
index 51c3aaa1772afff6e6d9ad580532a99485201861..f45588cfd84a9ed7288225811fbd4f6269ee1fc1 100644
--- a/FileIO/XmlIO/Qt/XMLQtInterface.cpp
+++ b/FileIO/XmlIO/Qt/XMLQtInterface.cpp
@@ -42,6 +42,7 @@ int XMLQtInterface::readFile(const QString &fileName)
 		return 0;
 	}
 	_fileData = file.readAll();
+	file.close();
 
 	if (!checkHash())
 		return 0;
@@ -110,7 +111,9 @@ bool XMLQtInterface::checkHash() const
 	QFile file(md5FileName);
 	if (file.open(QIODevice::ReadOnly))
 	{
-		if(file.readAll() == fileHash)
+		QByteArray referenceHash = file.readAll();
+		file.close();
+		if(referenceHash == fileHash)
 			return true;
 		INFO("Hashfile does not match data ... checking file ...");
 	}
diff --git a/FileIO/XmlIO/Qt/XmlLutReader.h b/FileIO/XmlIO/Qt/XmlLutReader.h
index 581412dbcd4d10f824e79288a6e0f7eb44929dd0..bc1ecbb9885ddf3c393c6ba46ad36f9868d6a40d 100644
--- a/FileIO/XmlIO/Qt/XmlLutReader.h
+++ b/FileIO/XmlIO/Qt/XmlLutReader.h
@@ -51,6 +51,7 @@ public:
 		if (docElement.nodeName().compare("ColorMap"))
 		{
 			ERR("XmlLutReader::readFromFile(): Unexpected XML root.");
+			file->close();
 			delete file;
 			return NULL;
 		}
@@ -97,6 +98,7 @@ public:
 
 		lut->SetTableRange(range[0], range[1]);
 
+		file->close();
 		delete file;
 
 		return lut;
diff --git a/GeoLib/Raster.cpp b/GeoLib/Raster.cpp
index 78168813a60a9094cad3811e83345414a338e041..c0623bd9c5bdb3ab33cc73f58f6879f117068b53 100644
--- a/GeoLib/Raster.cpp
+++ b/GeoLib/Raster.cpp
@@ -73,8 +73,8 @@ Raster* Raster::getRasterFromSurface(Surface const& sfc, double cell_size, doubl
 	Point const& ll (sfc.getAABB().getMinPoint());
 	Point const& ur (sfc.getAABB().getMaxPoint());
 
-	std::size_t n_cols = static_cast<size_t>(fabs(ur[0]-ll[0]) / cell_size)+1;
-	std::size_t n_rows = static_cast<size_t>(fabs(ur[1]-ll[1]) / cell_size)+1;
+	const std::size_t n_cols = static_cast<size_t>(fabs(ur[0]-ll[0]) / cell_size)+1;
+	const std::size_t n_rows = static_cast<size_t>(fabs(ur[1]-ll[1]) / cell_size)+1;
 	const size_t n_triangles (sfc.getNTriangles());
 	double *z_vals (new double[n_cols*n_rows]);
 	if (!z_vals) {
@@ -142,7 +142,7 @@ Raster* Raster::getRasterFromASCFile(std::string const& fname)
 		std::string s;
 		// read the data into the double-array
 		for (size_t j(0); j < n_rows; ++j) {
-			size_t idx ((n_rows - j - 1) * n_cols);
+			const size_t idx ((n_rows - j - 1) * n_cols);
 			for (size_t i(0); i < n_cols; ++i) {
 				in >> s;
 				values[idx+i] = strtod(BaseLib::replaceString(",", ".", s).c_str(),0);
@@ -215,18 +215,22 @@ Raster* Raster::getRasterFromSurferFile(std::string const& fname)
 
 	// header information
 	std::size_t n_cols(0), n_rows(0);
-	double xllcorner(0.0), yllcorner(0.0), cell_size(0.0), no_data_val(-9999);
+	double xllcorner(0.0), yllcorner(0.0), cell_size(0.0), min(0.0), max(0.0);
 
-	if (readSurferHeader(in, n_cols, n_rows, xllcorner, yllcorner, cell_size, no_data_val)) {
+	if (readSurferHeader(in, n_cols, n_rows, xllcorner, yllcorner, cell_size, min, max)) 
+	{
+		const double no_data_val (min-1);
 		double* values = new double[n_cols*n_rows];
 		std::string s;
 		// read the data into the double-array
-		for (size_t j(0); j < n_rows; ++j) {
-			size_t idx ((n_rows - j - 1) * n_cols);
-			for (size_t i(0); i < n_cols; ++i) {
+		for (size_t j(0); j < n_rows; ++j) 
+		{
+			const size_t idx (j * n_cols);
+			for (size_t i(0); i < n_cols; ++i) 
+			{
 				in >> s;
-				values[idx+i] = strtod(BaseLib::replaceString(",", ".", s).c_str(),0);
-
+				const double val (strtod(BaseLib::replaceString(",", ".", s).c_str(),0));
+				values[idx+i] = (val > max || val < min) ? no_data_val : val;
 			}
 		}
 		in.close();
@@ -241,11 +245,10 @@ Raster* Raster::getRasterFromSurferFile(std::string const& fname)
 }
 
 bool Raster::readSurferHeader(std::ifstream &in, size_t &n_cols, std::size_t &n_rows,
-				double &xllcorner, double &yllcorner, double &cell_size, double &no_data_val)
+				double &xllcorner, double &yllcorner, double &cell_size, double &min, double &max)
 {
 	std::string tag;
-	double min, max;
-
+	
 	in >> tag;
 
 	if (tag.compare("DSAA") != 0)
@@ -258,21 +261,19 @@ bool Raster::readSurferHeader(std::ifstream &in, size_t &n_cols, std::size_t &n_
 		in >> n_cols >> n_rows;
 		in >> min >> max;
 		xllcorner = min;
-		cell_size = (max-min)/(double)n_cols;
+		cell_size = (max-min)/static_cast<double>(n_cols);
 
 		in >> min >> max;
 		yllcorner = min;
 
-		if (ceil((max-min)/(double)n_rows) == ceil(cell_size))
+		if (ceil((max-min)/static_cast<double>(n_rows)) == ceil(cell_size))
 			cell_size = ceil(cell_size);
 		else
 		{
 			ERR("Error in readSurferHeader() - Anisotropic cellsize detected.");
 			return 0;
 		}
-		in >> min >> max; // ignore min- and max-values
-
-		no_data_val = 1.70141E+038;
+		in >> min >> max;
 	}
 
 	return true;
diff --git a/GeoLib/Raster.h b/GeoLib/Raster.h
index ff201534db62911fd0ab1b0ea975b7340d8b4a23..6e4c7acad9c30d5cb1a1c7964848f265b5616504 100644
--- a/GeoLib/Raster.h
+++ b/GeoLib/Raster.h
@@ -117,7 +117,7 @@ private:
 	 * \return True if the header could be read correctly, false otherwise.
 	 */
 	static bool readSurferHeader(std::ifstream &in, size_t &n_cols, std::size_t &n_rows,
-					double &xllcorner, double &yllcorner, double &cell_size, double &no_data_val);
+					double &xllcorner, double &yllcorner, double &cell_size, double &min, double &max);
 
 	void setCellSize(double cell_size);
 	void setNoDataVal (double no_data_val);
diff --git a/Gui/DataView/CMakeLists.txt b/Gui/DataView/CMakeLists.txt
index 661fc1fff5528694a73bd868a4028c2caa75b075..21a66b87756a4b3c9e3f9d0845b7a0ebbfc52dc6 100644
--- a/Gui/DataView/CMakeLists.txt
+++ b/Gui/DataView/CMakeLists.txt
@@ -28,7 +28,6 @@ set( SOURCES
 	ModelTreeItem.cpp
 	MeshLayerEditDialog.cpp
 	MshItem.cpp
-	MshLayerMapper.cpp
 	MshModel.cpp
 	MshQualitySelectionDialog.cpp
 	MshTabWidget.cpp
@@ -98,7 +97,6 @@ set( HEADERS
 	GeoTreeItem.h
 	ModelTreeItem.h
 	MshItem.h
-	MshLayerMapper.h
 	ProcessItem.h
 )
 
diff --git a/Gui/DataView/CondFromRasterDialog.cpp b/Gui/DataView/CondFromRasterDialog.cpp
index 538e1b3c4a097a1f5eb63be2893040324f18f781..5c1da975ae053353c3ad963571f73f1c1d0eae13 100644
--- a/Gui/DataView/CondFromRasterDialog.cpp
+++ b/Gui/DataView/CondFromRasterDialog.cpp
@@ -52,15 +52,15 @@ void CondFromRasterDialog::on_selectButton_pressed()
 	QString geotiffExtension("");
 #endif
 	QString fileName = QFileDialog::getOpenFileName(this, "Select raster file",
-					settings.value("lastOpenedConditionsFileDirectory").toString(), QString(
-									"Raster files (*.asc *.grd);;") .arg(geotiffExtension));
+	                                                settings.value("lastOpenedRasterFileDirectory").toString(), 
+	                                                QString("Raster files (*.asc *.grd);;").arg(geotiffExtension));
 
 	if (!fileName.isEmpty())
 	{
 		this->rasterEdit->setText(fileName);
 
-		QDir dir = QDir(fileName);
-		settings.setValue("lastOpenedConditionsFileDirectory", dir.absolutePath());
+		QFileInfo fi(fileName);
+		settings.setValue("lastOpenedRasterFileDirectory", fi.absolutePath());
 	}
 }
 
diff --git a/Gui/DataView/MeshElementRemovalDialog.cpp b/Gui/DataView/MeshElementRemovalDialog.cpp
index 11397168b05f9867c049c531c2464e1ce3073187..0c855565cc826089bf148bf5bbdd3501161a841b 100644
--- a/Gui/DataView/MeshElementRemovalDialog.cpp
+++ b/Gui/DataView/MeshElementRemovalDialog.cpp
@@ -100,7 +100,7 @@ void MeshElementRemovalDialog::accept()
 			emit meshAdded(new_mesh);
 		else
 		{
-			unsigned error_code (ex.getErrorCode());
+			const unsigned error_code (ex.getErrorCode());
 			if (error_code == 1)
 				OGSError::box("The current selection removes ALL mesh elements.\nPlease change the selection.");
 			if (error_code == 2)
diff --git a/Gui/DataView/MeshLayerEditDialog.cpp b/Gui/DataView/MeshLayerEditDialog.cpp
index 9db1c4b653222f6239492fc371493e418f5b60e6..5692a07ce325ae833b948273d37df2874fcb4cab 100644
--- a/Gui/DataView/MeshLayerEditDialog.cpp
+++ b/Gui/DataView/MeshLayerEditDialog.cpp
@@ -15,11 +15,15 @@
 
 #include "MeshLayerEditDialog.h"
 
+// ThirdParty/logog
+#include "logog/include/logog.hpp"
+
 #include "OGSError.h"
 #include "StringTools.h"
 #include "Mesh.h"
 
 #include <QCheckBox>
+#include <QFileInfo>
 #include <QFileDialog>
 #include <QPushButton>
 #include <QSettings>
@@ -188,7 +192,7 @@ void MeshLayerEditDialog::accept()
 				const std::string imgPath ( this->_edits[0]->text().toStdString() );
 				const double noDataReplacementValue = this->_noDataReplacementEdit->text().toDouble();
 				if (!imgPath.empty())
-					result = MshLayerMapper::LayerMapping(new_mesh, imgPath, nLayers, 0, noDataReplacementValue);
+					result = MeshLayerMapper::LayerMapping(*new_mesh, imgPath, nLayers, 0, noDataReplacementValue);
 			}
 			else
 			{
@@ -197,12 +201,11 @@ void MeshLayerEditDialog::accept()
 				{
 					// "100" is just a default size to have any value for extruding 2D elements.
 					// The actual mapping based on a raster file will be performed later.
-					float thickness = (_use_rasters) ? 100 : (this->_edits[i]->text().toFloat());
-					if (thickness > std::numeric_limits<float>::epsilon())
-						layer_thickness.push_back(thickness);
+					const float thickness = (_use_rasters) ? 100 : (this->_edits[i]->text().toFloat());
+					layer_thickness.push_back(thickness);
 				}
 
-				new_mesh = MshLayerMapper::CreateLayers(_msh, layer_thickness);
+				new_mesh = MeshLayerMapper::CreateLayers(*_msh, layer_thickness);
 
 				if (_use_rasters)
 				{
@@ -212,13 +215,13 @@ void MeshLayerEditDialog::accept()
 						const double noDataReplacement = (i==0) ? 0.0 : -9999.0;
 						if (!imgPath.empty())
 						{
-							result = MshLayerMapper::LayerMapping(new_mesh, imgPath, nLayers, i, noDataReplacement);
+							result = MeshLayerMapper::LayerMapping(*new_mesh, imgPath, nLayers, i, noDataReplacement);
 							if (result==0) break;
 						}
 					}
 					if (this->_edits[0]->text().length()>0)
 					{
-						MeshLib::Mesh* final_mesh = MshLayerMapper::blendLayersWithSurface(new_mesh, nLayers, this->_edits[0]->text().toStdString());
+						MeshLib::Mesh* final_mesh = MeshLayerMapper::blendLayersWithSurface(*new_mesh, nLayers, this->_edits[0]->text().toStdString());
 						delete new_mesh;
 						new_mesh = final_mesh;
 					}
@@ -250,9 +253,9 @@ void MeshLayerEditDialog::getFileName()
 	QPushButton* button = dynamic_cast<QPushButton*>(this->sender());
 	QSettings settings;
 	QString filename = QFileDialog::getOpenFileName(this, "Select raster file to open",
-	                                                settings.value("lastOpenedFileDirectory").toString(),
+	                                                settings.value("lastOpenedRasterFileDirectory").toString(),
 	                                                "ASCII raster files (*.asc);;All files (* *.*)");
 	_fileButtonMap[button]->setText(filename);
-	QDir dir = QDir(filename);
-	settings.setValue("lastOpenedFileDirectory", dir.absolutePath());
+	QFileInfo fi(filename);
+	settings.setValue("lastOpenedRasterFileDirectory", fi.absolutePath());
 }
diff --git a/Gui/DataView/MeshLayerEditDialog.h b/Gui/DataView/MeshLayerEditDialog.h
index 1b02632fb182eced09000282a8ec44f8d92e767e..89fabac92366e00ed6a5e243e4f2f92d362ffa1c 100644
--- a/Gui/DataView/MeshLayerEditDialog.h
+++ b/Gui/DataView/MeshLayerEditDialog.h
@@ -18,7 +18,7 @@
 #include "ui_MeshLayerEdit.h"
 #include <QDialog>
 
-#include "MshLayerMapper.h"
+#include "MeshGenerators/MeshLayerMapper.h"
 
 class QPushButton;
 class QCheckBox;
diff --git a/Gui/VtkVis/VtkCompositeColormapToImageFilter.cpp b/Gui/VtkVis/VtkCompositeColormapToImageFilter.cpp
index fb0abb34f45b466447ba9523f64f2a2225b5a19f..f9224336eeb8cd2bee095065d6d1c1f429bfa252 100644
--- a/Gui/VtkVis/VtkCompositeColormapToImageFilter.cpp
+++ b/Gui/VtkVis/VtkCompositeColormapToImageFilter.cpp
@@ -52,7 +52,7 @@ void VtkCompositeColormapToImageFilter::init()
 	double range[2];
 	dynamic_cast<vtkImageAlgorithm*>(_inputAlgorithm)->GetOutput()->GetPointData()->GetScalars()->GetRange(range);
 
-	if (!fileName.length()==0)
+	if (fileName.length() > 0)
 	{
 		colormap = FileIO::XmlLutReader::readFromFile(fileName);
 		settings.setValue("lastOpenedLookupTableFileDirectory", fileName);
diff --git a/Gui/VtkVis/VtkRaster.cpp b/Gui/VtkVis/VtkRaster.cpp
index b31902955fd54762de786319aa30ed66a21683b4..fc50caa9b9522c9c044df19bb59e3a23f91686c4 100644
--- a/Gui/VtkVis/VtkRaster.cpp
+++ b/Gui/VtkVis/VtkRaster.cpp
@@ -14,6 +14,7 @@
 
 #include "VtkRaster.h"
 
+#include <algorithm>
 #include <cmath>
 #include <iomanip>
 #include <iostream>
@@ -85,24 +86,17 @@ vtkImageImport* VtkRaster::loadImageFromArray(double const*const data_array, dou
 {
 	const unsigned length = height*width;
 	float* data = new float[length*2];
-	float max_val=noData;
+	float max_val = *std::max_element(data_array, data_array+length);
 	for (unsigned j=0; j<length; ++j)
 	{
 		data[j*2] = static_cast<float>(data_array[j]);
-		max_val = (data[j*2]>max_val) ? data[j*2] : max_val;
-	}
-	for (unsigned j=0; j<length; ++j)
-	{
-		if (data[j*2]==noData)
+		if (fabs(data[j*2]-noData) < std::numeric_limits<double>::epsilon())
 		{
 			data[j*2] = max_val;
 			data[j*2+1] = 0;
 		}
 		else
-		{
-			//data[j*2] = max_val-data[j*2];//delete;
 			data[j*2+1] = max_val;
-		}
 	}
 
 	vtkImageImport* image = vtkImageImport::New();
@@ -119,9 +113,6 @@ vtkImageImport* VtkRaster::loadImageFromArray(double const*const data_array, dou
 	return image;
 }
 
-
-
-
 #ifdef GEOTIFF_FOUND
 vtkImageImport* VtkRaster::loadImageFromTIFF(const std::string &fileName,
                                   double &x0, double &y0,
diff --git a/Gui/mainwindow.cpp b/Gui/mainwindow.cpp
index 8bef3704079063261b0484de3bacf7529735f9cb..33a9dba8cd96cf5aeb3e2ca7be998759ca614644 100644
--- a/Gui/mainwindow.cpp
+++ b/Gui/mainwindow.cpp
@@ -514,7 +514,7 @@ void MainWindow::save()
 void MainWindow::loadFile(ImportFileType::type t, const QString &fileName)
 {
 	QFile file(fileName);
-	if (!file.open(QFile::ReadOnly))
+	if (!file.exists())
 	{
 		QMessageBox::warning(this, tr("Application"), tr(
 		                             "Cannot read file %1:\n%2.").arg(fileName).arg(
diff --git a/MathLib/MathTools.cpp b/MathLib/MathTools.cpp
index 7de70163cbea2329342798b05848b02c8f505e3c..651fa75369b10a4218753378150339db1694d346 100644
--- a/MathLib/MathTools.cpp
+++ b/MathLib/MathTools.cpp
@@ -86,13 +86,4 @@ double calcTetrahedronVolume(const double* x1, const double* x2, const double* x
 	          + (x1[2] - x4[2]) * ((x2[0] - x4[0]) * (x3[1] - x4[1]) - (x2[1] - x4[1]) * (x3[0] - x4[0]))) / 6.0;
 }
 
-void MPhi2D(double* vf, double r, double s)
-{
-	vf[0] = (1.0 + r) * (1.0 + s);
-	vf[1] = (1.0 - r) * (1.0 + s);
-	vf[2] = (1.0 - r) * (1.0 - s);
-	vf[3] = (1.0 + r) * (1.0 - s);
-	for (unsigned i = 0; i < 4; i++)
-		vf[i] *= 0.25;
-}
 } // namespace
diff --git a/MathLib/MathTools.h b/MathLib/MathTools.h
index f205d25b598c700d82d2edbdf60929117d9c3f22..585d141b14c405b6143eeb70d2ac50a59fc559ab 100644
--- a/MathLib/MathTools.h
+++ b/MathLib/MathTools.h
@@ -177,8 +177,6 @@ T fastpow (T base, std::size_t exp)
 	return result;
 }
 
-/// 2D linear interpolation function (TODO adopted from geo_mathlib)
-void MPhi2D(double* vf, double r, double s);
 } // namespace
 
 #endif /* MATHTOOLS_H_ */
diff --git a/MeshLib/MeshEditing/ElementExtraction.cpp b/MeshLib/MeshEditing/ElementExtraction.cpp
index cbdd542d9bb26ec563906bffa5572a37ba401dbc..19e31c305c38f10e98c2ffa69771dfe18223b6f0 100644
--- a/MeshLib/MeshEditing/ElementExtraction.cpp
+++ b/MeshLib/MeshEditing/ElementExtraction.cpp
@@ -15,6 +15,7 @@
 #include "ElementExtraction.h"
 #include "Mesh.h"
 #include "Elements/Element.h"
+#include "MeshEditing/DuplicateMeshComponents.h"
 #include "AABB.h"
 
 #include "logog/include/logog.hpp"
@@ -42,9 +43,8 @@ MeshLib::Mesh* ElementExtraction::removeMeshElements(const std::string &new_mesh
 	INFO("Removing total %d elements...", _marked_elements.size());
 	std::vector<MeshLib::Element*> tmp_elems = excludeElements(_mesh.getElements(), _marked_elements);
 	INFO("%d elements remain in mesh.", tmp_elems.size());
-	std::vector<MeshLib::Node*> new_nodes;
-	std::vector<MeshLib::Element*> new_elems;
-	copyNodesElements(_mesh.getNodes(), tmp_elems, new_nodes, new_elems);
+	std::vector<MeshLib::Node*> new_nodes = MeshLib::copyNodeVector(_mesh.getNodes());
+	std::vector<MeshLib::Element*> new_elems = MeshLib::copyElementVector(tmp_elems, new_nodes);
 
 	// create a new mesh object. Unsued nodes are removed during construction
 	if (!new_elems.empty())
@@ -57,7 +57,7 @@ MeshLib::Mesh* ElementExtraction::removeMeshElements(const std::string &new_mesh
 	}
 }
 
-void ElementExtraction::searchByMaterialID(unsigned matID)
+std::size_t ElementExtraction::searchByMaterialID(unsigned matID)
 {
 	const std::vector<MeshLib::Element*> &ele_vec (this->_mesh.getElements());
 	std::vector<std::size_t> matchedIDs;
@@ -68,9 +68,10 @@ void ElementExtraction::searchByMaterialID(unsigned matID)
 		i++;
 	}
 	this->updateUnion(matchedIDs);
+	return matchedIDs.size();
 }
 
-void ElementExtraction::searchByElementType(MeshElemType eleType)
+std::size_t ElementExtraction::searchByElementType(MeshElemType eleType)
 {
 	const std::vector<MeshLib::Element*> &ele_vec (this->_mesh.getElements());
 	std::vector<std::size_t> matchedIDs;
@@ -81,9 +82,10 @@ void ElementExtraction::searchByElementType(MeshElemType eleType)
 		i++;
 	}
 	this->updateUnion(matchedIDs);
+	return matchedIDs.size();
 }
 
-void ElementExtraction::searchByZeroContent()
+std::size_t ElementExtraction::searchByZeroContent()
 {
 	const std::vector<MeshLib::Element*> &ele_vec (this->_mesh.getElements());
 	std::vector<std::size_t> matchedIDs;
@@ -94,9 +96,10 @@ void ElementExtraction::searchByZeroContent()
 		i++;
 	}
 	this->updateUnion(matchedIDs);
+	return matchedIDs.size();
 }
 
-void ElementExtraction::searchByBoundingBox(const MeshLib::Node &x1, const MeshLib::Node &x2)
+std::size_t ElementExtraction::searchByBoundingBox(const MeshLib::Node &x1, const MeshLib::Node &x2)
 {
 	const std::vector<MeshLib::Element*> &ele_vec (this->_mesh.getElements());
 	std::vector<MeshLib::Node> extent;
@@ -117,6 +120,7 @@ void ElementExtraction::searchByBoundingBox(const MeshLib::Node &x1, const MeshL
 			}
 	}
 	this->updateUnion(matchedIDs);
+	return matchedIDs.size();
 }
 
 void ElementExtraction::updateUnion(const std::vector<std::size_t> &vec)
@@ -143,29 +147,4 @@ std::vector<MeshLib::Element*> ElementExtraction::excludeElements(const std::vec
 	return vec_dest_eles;
 }
 
-void ElementExtraction::copyNodesElements(	const std::vector<MeshLib::Node*> &src_nodes, const std::vector<MeshLib::Element*> &src_elems,
-						std::vector<MeshLib::Node*> &dst_nodes, std::vector<MeshLib::Element*> &dst_elems) const 
-{
-	// 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_elems.resize(src_elems.size());
-	for (std::size_t i=0; i<dst_elems.size(); i++) {
-		auto* src_elem = src_elems[i];
-		auto* dst_elem = src_elem->clone();
-		for (unsigned j=0; j<src_elem->getNNodes(); j++) {
-			dst_elem->setNode(j, dst_nodes[src_elem->getNode(j)->getID()]);
-		}
-		dst_elems[i] = dst_elem;
-	}
-}
-
-
-
-
-
 } // end namespace MeshLib
diff --git a/MeshLib/MeshEditing/ElementExtraction.h b/MeshLib/MeshEditing/ElementExtraction.h
index eb310a9504e4f861ac5eac2fcb2a556d92b9cdb7..eeb9b04e5e3e1b2c9b8aa82db0be756ce07a9b7d 100644
--- a/MeshLib/MeshEditing/ElementExtraction.h
+++ b/MeshLib/MeshEditing/ElementExtraction.h
@@ -41,16 +41,16 @@ public:
 	MeshLib::Mesh* removeMeshElements(const std::string &new_mesh_name);
 	
 	/// Marks all elements with the given Material ID.
-	void searchByMaterialID(unsigned matID);
+	std::size_t searchByMaterialID(unsigned matID);
 
 	/// Marks all elements of the given element type.
-	void searchByElementType(MeshElemType eleType);
+	std::size_t searchByElementType(MeshElemType eleType);
 
 	/// Marks all elements with a volume smaller than std::numeric_limits<double>::epsilon().
-	void searchByZeroContent();
+	std::size_t searchByZeroContent();
 
 	/// Marks all elements with at least one node outside the bounding box spanned by x1 and x2;
-	void searchByBoundingBox(const MeshLib::Node &x1, const MeshLib::Node &x2);
+	std::size_t searchByBoundingBox(const MeshLib::Node &x1, const MeshLib::Node &x2);
 	
 
 private:
@@ -60,10 +60,6 @@ private:
 	/// Removes elements from vec_removed in vec_src_elems
 	std::vector<MeshLib::Element*> excludeElements(const std::vector<MeshLib::Element*> & vec_src_elems, const std::vector<std::size_t> &vec_removed) const;
 
-	/// Copies nodes and elements of the original mesh for constructing the new mesh
-	void copyNodesElements(const std::vector<MeshLib::Node*> &src_nodes, const std::vector<MeshLib::Element*> &src_elems,
-						   std::vector<MeshLib::Node*> &dst_nodes, std::vector<MeshLib::Element*> &dst_elems) const;
-
 	/// The mesh from which elements should be removed.
 	const MeshLib::Mesh &_mesh;
 	/// The vector of element indices that should be removed.
diff --git a/MeshLib/MeshEditing/MeshRevision.cpp b/MeshLib/MeshEditing/MeshRevision.cpp
index 3a120354e627619538bbb403d2725ebc7efb1971..3f9de650f1805234b216ff166dab406db8055ea7 100644
--- a/MeshLib/MeshEditing/MeshRevision.cpp
+++ b/MeshLib/MeshEditing/MeshRevision.cpp
@@ -45,12 +45,8 @@ MeshRevision::MeshRevision(MeshLib::Mesh &mesh) :
 
 MeshLib::Mesh* MeshRevision::collapseNodes(const std::string &new_mesh_name, double eps)
 {
-	std::vector<MeshLib::Node*> new_nodes = this->constructNewNodesArray(this->collapseNodeIndeces(eps));
-	std::vector<MeshLib::Element*> new_elements;
-	new_elements.reserve(this->_mesh.getNElements());
-	std::vector<MeshLib::Element*> const& elements(this->_mesh.getElements());
-	for (auto elem = elements.begin(); elem != elements.end(); ++elem)
-		new_elements.push_back(MeshLib::copyElement(*elem, new_nodes));
+	std::vector<MeshLib::Node*> new_nodes (this->constructNewNodesArray(this->collapseNodeIndeces(eps)));
+	std::vector<MeshLib::Element*> new_elements (MeshLib::copyElementVector(_mesh.getElements(), new_nodes));
 	this->resetNodeIDs();
 	return new MeshLib::Mesh(new_mesh_name, new_nodes, new_elements);
 }
diff --git a/MeshLib/MeshEditing/MeshRevision.h b/MeshLib/MeshEditing/MeshRevision.h
index 1960b69fd481cb36d2f676f81678e3c94f4904ce..3e1edacb4450a4c913cfe32b3f0927b54306134f 100644
--- a/MeshLib/MeshEditing/MeshRevision.h
+++ b/MeshLib/MeshEditing/MeshRevision.h
@@ -35,6 +35,7 @@ namespace MeshLib {
  * reduces elements accordingly.
  */
 class MeshRevision {
+
 public:
 	/** 
 	 * Constructor
diff --git a/MeshLib/MeshGenerators/MeshGenerator.cpp b/MeshLib/MeshGenerators/MeshGenerator.cpp
index 5a9e51976f830b32203ec2e62ebfca002ced36f7..cd6576856bab399266c7d8e1183e391d04aea06c 100644
--- a/MeshLib/MeshGenerators/MeshGenerator.cpp
+++ b/MeshLib/MeshGenerators/MeshGenerator.cpp
@@ -68,10 +68,10 @@ Mesh* MeshGenerator::generateRegularQuadMesh(
 }
 
 Mesh* MeshGenerator::generateRegularQuadMesh(const unsigned n_x_cells,
-	                          const unsigned n_y_cells,
-                              const double   cell_size,
+                              const unsigned n_y_cells,
+                              const double cell_size,
                               GeoLib::Point const& origin,
-							  std::string   const& mesh_name)
+                              std::string const& mesh_name)
 {
 	//nodes
 	const unsigned n_x_nodes (n_x_cells+1);
@@ -79,7 +79,7 @@ Mesh* MeshGenerator::generateRegularQuadMesh(const unsigned n_x_cells,
 	std::vector<Node*> nodes;
 	nodes.reserve(n_x_nodes * n_y_nodes);	
 
-	for (std::size_t i = 0, node_id = 0; i < n_y_nodes; i++)
+	for (std::size_t i = 0; i < n_y_nodes; i++)
 	{
 		const double y_offset (origin[1] + cell_size * i);
 		for (std::size_t j = 0; j < n_x_nodes; j++)
@@ -130,7 +130,7 @@ Mesh* MeshGenerator::generateRegularHexMesh(const unsigned n_x_cells,
 	std::vector<Node*> nodes;
 	nodes.reserve(n_x_nodes * n_y_nodes * n_z_nodes);
 
-	for (std::size_t i = 0, node_id = 0; i < n_z_nodes; i++)
+	for (std::size_t i = 0; i < n_z_nodes; i++)
 	{
 		const double z_offset (origin[2] + cell_size * i);
 		for (std::size_t j = 0; j < n_y_nodes; j++)
diff --git a/Gui/DataView/MshLayerMapper.cpp b/MeshLib/MeshGenerators/MeshLayerMapper.cpp
similarity index 50%
rename from Gui/DataView/MshLayerMapper.cpp
rename to MeshLib/MeshGenerators/MeshLayerMapper.cpp
index 24b3731ffdabcfe21321966ed758c29973f0a686..59c600ef2878f7bcd62adf4b4e1e5815003a382c 100644
--- a/Gui/DataView/MshLayerMapper.cpp
+++ b/MeshLib/MeshGenerators/MeshLayerMapper.cpp
@@ -1,8 +1,8 @@
 /**
- * \file
+ * \file   MeshLayerMapper.cpp
  * \author Karsten Rink
  * \date   2010-11-01
- * \brief  Implementation of the MshLayerMapper class.
+ * \brief  Implementation of the MeshLayerMapper class.
  *
  * \copyright
  * Copyright (c) 2013, OpenGeoSys Community (http://www.opengeosys.org)
@@ -18,7 +18,7 @@
 // ThirdParty/logog
 #include "logog/include/logog.hpp"
 
-#include "MshLayerMapper.h"
+#include "MeshLayerMapper.h"
 // GeoLib
 #include "Raster.h"
 
@@ -33,199 +33,179 @@
 #include "MathTools.h"
 
 
-MeshLib::Mesh* MshLayerMapper::CreateLayers(const MeshLib::Mesh* mesh, const std::vector<float> &thickness)
+MeshLib::Mesh* MeshLayerMapper::CreateLayers(const MeshLib::Mesh &mesh, const std::vector<float> &layer_thickness_vector)
 {
-	std::size_t nLayers(thickness.size());
-	bool throw_error(false);
-	for (unsigned i=0; i<nLayers; ++i)
-		if (thickness[i]<=0)
-			throw_error = true;
-	if (nLayers < 1 || mesh->getDimension() != 2)
-		throw_error = true;
-
-	if (throw_error)
+	std::vector<float> thickness;
+	for (std::size_t i=0; i<layer_thickness_vector.size(); ++i)
+		if (layer_thickness_vector[i] > std::numeric_limits<float>::epsilon())
+			thickness.push_back(layer_thickness_vector[i]);
+		else
+			WARN ("Ignoring layer %d with thickness %f.", i, layer_thickness_vector[i]);
+
+	const std::size_t nLayers(thickness.size());
+	if (nLayers < 1 || mesh.getDimension() != 2)
 	{
 		ERR("MshLayerMapper::CreateLayers(): A 2D mesh with nLayers > 0 is required as input.");
 		return nullptr;
 	}
 
-	const size_t nNodes = mesh->getNNodes();
+	const size_t nNodes = mesh.getNNodes();
 	// count number of 2d elements in the original mesh
-	const std::size_t nElems (std::count_if(mesh->getElements().begin(), mesh->getElements().end(),
+	const std::size_t nElems (std::count_if(mesh.getElements().begin(), mesh.getElements().end(),
 			[](MeshLib::Element const* elem) { return (elem->getDimension() == 2);}));
 
-	const std::vector<MeshLib::Node*> nodes = mesh->getNodes();
-	const std::vector<MeshLib::Element*> elems = mesh->getElements();
-	std::vector<MeshLib::Node*> new_nodes(nNodes + (nLayers * nNodes));
-	std::vector<MeshLib::Element*> new_elems(nElems * nLayers);
-	double z_offset(0.0);
+	const std::size_t nOrgElems (mesh.getNElements());
+	const std::vector<MeshLib::Node*> &nodes = mesh.getNodes();
+	const std::vector<MeshLib::Element*> &elems = mesh.getElements();
+	std::vector<MeshLib::Node*> new_nodes;
+	new_nodes.reserve(nNodes + (nLayers * nNodes));
+	std::vector<MeshLib::Element*> new_elems;
+	new_elems.reserve(nElems * nLayers);
+	double z_offset (0.0);
 
 	for (unsigned layer_id = 0; layer_id <= nLayers; ++layer_id)
 	{
 		// add nodes for new layer
 		unsigned node_offset (nNodes * layer_id);
-		unsigned elem_offset (nElems * (layer_id-1));
-		if (layer_id>0) z_offset += thickness[layer_id-1];
+		if (layer_id > 0) z_offset += thickness[layer_id-1];
 		for (unsigned i = 0; i < nNodes; ++i)
 		{
 			const double* coords = nodes[i]->getCoords();
-			new_nodes[node_offset+i] = new MeshLib::Node(coords[0], coords[1], coords[2]-z_offset, node_offset+i);
+			new_nodes.push_back (new MeshLib::Node(coords[0], coords[1], coords[2]-z_offset, node_offset+i));
 		}
 
 		// starting with 2nd layer create prism or hex elements connecting the last layer with the current one
 		if (layer_id > 0)
 		{
-			if (thickness[layer_id-1] > 0)
+			node_offset -= nNodes;
+			const unsigned mat_id (nLayers - layer_id);
+
+			for (unsigned i = 0; i < nOrgElems; ++i)
 			{
-				node_offset -= nNodes;
-				const unsigned mat_id (nLayers - layer_id);
+				const MeshLib::Element* sfc_elem( elems[i] );
+				if (sfc_elem->getDimension() != 2)
+				{
+					WARN("MshLayerMapper::CreateLayers() - Method can only handle 2D mesh elements.");
+					WARN("Skipping Element %d of type \"%s\".", i, MeshElemType2String(sfc_elem->getGeomType()).c_str());
+					continue;
+				}
+				
+				const unsigned nElemNodes(sfc_elem->getNNodes());
+				MeshLib::Node** e_nodes = new MeshLib::Node*[2*nElemNodes];
 
-				// counts the 2d elements (within the layer),
-				// used as a part of index computation in new_elems
-				std::size_t cnt(0);
-				for (unsigned i = 0; i < mesh->getNElements(); ++i)
+				for (unsigned j=0; j<nElemNodes; ++j)
 				{
-					const MeshLib::Element* sfc_elem( elems[i] );
-					if (sfc_elem->getDimension() == 2)
-					{
-						const unsigned nElemNodes(sfc_elem->getNNodes());
-						MeshLib::Node** e_nodes = new MeshLib::Node*[2*nElemNodes];
-
-						for (unsigned j=0; j<nElemNodes; ++j)
-						{
-							const unsigned node_id = sfc_elem->getNode(j)->getID() + node_offset;
-							e_nodes[j] = new_nodes[node_id+nNodes];
-							e_nodes[j+nElemNodes] = new_nodes[node_id];
-						}
-						if (sfc_elem->getGeomType() == MeshElemType::TRIANGLE)	// extrude triangles to prism
-							new_elems[elem_offset+cnt] = new MeshLib::Prism(e_nodes, mat_id);
-						else if (sfc_elem->getGeomType() == MeshElemType::QUAD)	// extrude quads to hexes
-							new_elems[elem_offset+cnt] = new MeshLib::Hex(e_nodes, mat_id);
-						cnt++;
-					}
-					else
-					{
-						WARN("MshLayerMapper::CreateLayers() - Method can only handle 2D mesh elements.");
-						WARN("Skipping Element %d of type \"%s\".", i, MeshElemType2String(sfc_elem->getGeomType()).c_str());
-					}
+					const unsigned node_id = sfc_elem->getNode(j)->getID() + node_offset;
+					e_nodes[j] = new_nodes[node_id+nNodes];
+					e_nodes[j+nElemNodes] = new_nodes[node_id];
 				}
-			}
-			else
-			{
-				ERR("Error in MshLayerMapper::CreateLayers() - Layer thickness for layer %d is %f (needs to be >0).", (layer_id-1), thickness[layer_id-1]);
-				return nullptr;
+				if (sfc_elem->getGeomType() == MeshElemType::TRIANGLE)	// extrude triangles to prism
+					new_elems.push_back (new MeshLib::Prism(e_nodes, mat_id));
+				else if (sfc_elem->getGeomType() == MeshElemType::QUAD)	// extrude quads to hexes
+					new_elems.push_back (new MeshLib::Hex(e_nodes, mat_id));
 			}
 		}
 	}
 	return new MeshLib::Mesh("SubsurfaceMesh", new_nodes, new_elems);
 }
 
-int MshLayerMapper::LayerMapping(MeshLib::Mesh* new_mesh, const std::string &rasterfile,
+int MeshLayerMapper::LayerMapping(MeshLib::Mesh &new_mesh, const std::string &rasterfile,
                                  const unsigned nLayers, const unsigned layer_id, double noDataReplacementValue = 0.0)
 {
-	if (new_mesh == nullptr)
+	if (nLayers < layer_id)
 	{
-		ERR("MshLayerMapper::LayerMapping() - Passed Mesh is NULL.");
+		ERR("MshLayerMapper::LayerMapping() - Mesh has only %d Layers, cannot assign layer %d.", nLayers, layer_id);
 		return 0;
 	}
 
-	if (nLayers >= layer_id)
+	const GeoLib::Raster *raster(GeoLib::Raster::getRasterFromASCFile(rasterfile));
+	if (! raster) {
+		ERR("MshLayerMapper::LayerMapping - could not read raster file %s", rasterfile.c_str());
+		return 0;
+	}
+	const double x0(raster->getOrigin()[0]);
+	const double y0(raster->getOrigin()[1]);
+	const double delta(raster->getRasterPixelDistance());
+	const double no_data(raster->getNoDataValue());
+	const std::size_t width(raster->getNCols());
+	const std::size_t height(raster->getNRows());
+	double const*const elevation(raster->begin());
+
+	const std::pair<double, double> xDim(x0, x0 + width * delta); // extension in x-dimension
+	const std::pair<double, double> yDim(y0, y0 + height * delta); // extension in y-dimension
+
+	const size_t nNodes (new_mesh.getNNodes());
+	const size_t nNodesPerLayer (nNodes / (nLayers+1));
+
+	const size_t firstNode (layer_id * nNodesPerLayer);
+	const size_t lastNode  (firstNode + nNodesPerLayer);
+
+	const double half_delta = 0.5*delta;
+	const std::vector<MeshLib::Node*> &nodes = new_mesh.getNodes();
+	for (unsigned i = firstNode; i < lastNode; ++i)
 	{
-		GeoLib::Raster *raster(GeoLib::Raster::getRasterFromASCFile(rasterfile));
-		if (! raster) {
-			ERR("MshLayerMapper::LayerMapping - could not read raster file %s", rasterfile.c_str());
-			return 0;
-		}
-		const double x0(raster->getOrigin()[0]);
-		const double y0(raster->getOrigin()[1]);
-		const double delta(raster->getRasterPixelDistance());
-		const double no_data(raster->getNoDataValue());
-		const std::size_t width(raster->getNCols());
-		const std::size_t height(raster->getNRows());
-		double const*const elevation(raster->begin());
-
-		const std::pair<double, double> xDim(x0, x0 + width * delta); // extension in x-dimension
-		const std::pair<double, double> yDim(y0, y0 + height * delta); // extension in y-dimension
+		const double* coords (nodes[i]->getCoords());
 
-		const size_t nNodes = new_mesh->getNNodes();
-		const size_t nNodesPerLayer = nNodes / (nLayers+1);
-
-		const size_t firstNode = layer_id * nNodesPerLayer;
-		const size_t lastNode  = firstNode + nNodesPerLayer;
-
-		const double half_delta = 0.5*delta;
-		const std::vector<MeshLib::Node*> nodes = new_mesh->getNodes();
-		for (unsigned i = firstNode; i < lastNode; ++i)
+		if (!isNodeOnRaster(*nodes[i], xDim, yDim))
 		{
-			const double* coords = nodes[i]->getCoords();
+			// use either default value or elevation from layer above
+			nodes[i]->updateCoordinates(coords[0], coords[1], noDataReplacementValue);
+			continue;
+		}
 
-			if (!isNodeOnRaster(*nodes[i], xDim, yDim))
+		// position in raster
+		const double xPos ((coords[0] - xDim.first) / delta);
+		const double yPos ((coords[1] - yDim.first) / delta);
+		// raster cell index
+		const size_t xIdx (static_cast<size_t>(floor(xPos)));
+		const size_t yIdx (static_cast<size_t>(floor(yPos)));
+
+		// weights for bilinear interpolation		
+		const double xShift = fabs(xPos-(xIdx+half_delta))/delta;
+		const double yShift = fabs(yPos-(yIdx+half_delta))/delta;
+		std::array<double,4> weight = { (1-xShift)*(1-xShift), xShift*(1-yShift), xShift*yShift, (1-xShift)*yShift };
+
+		// neightbors to include in interpolation
+		const int xShiftIdx = (xPos-xIdx-half_delta>=0) ? 1 : -1;
+		const int yShiftIdx = (yPos-yIdx-half_delta>=0) ? 1 : -1;
+		const std::array<int,4> x_nb = { 0, xShiftIdx, xShiftIdx, 0 };
+		const std::array<int,4> y_nb = { 0, 0, yShiftIdx, yShiftIdx };
+
+		// get pixel values
+		std::array<double,4>  pix_val;
+		unsigned no_data_count (0);
+		for (unsigned j=0; j<4; ++j)
+		{
+			pix_val[j] = elevation[(yIdx + y_nb[j])*width + (xIdx + x_nb[j])];
+			if (fabs(pix_val[j] - no_data) < std::numeric_limits<double>::epsilon())
 			{
-				if (layer_id == 0) // use default value
-					nodes[i]->updateCoordinates(coords[0], coords[1], noDataReplacementValue);
-				else // use z-value from layer above
-					nodes[i]->updateCoordinates(coords[0], coords[1], (*nodes[i-nNodesPerLayer])[2]);
-				continue;
+				weight[j] = 0;
+				no_data_count++;
 			}
+		}
 
-			// position in raster
-			const double xPos ((coords[0] - xDim.first) / delta);
-			const double yPos ((coords[1] - yDim.first) / delta);
-			// raster cell index
-			const size_t xIdx (static_cast<size_t>(floor(xPos)));
-			const size_t yIdx (static_cast<size_t>(floor(yPos)));
-
-			// deviation of mesh node from centre of raster cell ( in [-1:1) because it is normalised by delta/2 )
-			const double xShift = (xPos-xIdx-half_delta)/half_delta;
-			const double yShift = (yPos-yIdx-half_delta)/half_delta;
-
-			const int xShiftIdx = static_cast<int>((xShift>=0) ? ceil(xShift) : floor(xShift));
-			const int yShiftIdx = static_cast<int>((yShift>=0) ? ceil(yShift) : floor(yShift));
-
-			// determining the neighbouring pixels that add weight to the interpolation
-			const int x_nb[4] = {0, xShiftIdx, xShiftIdx, 0};
-			const int y_nb[4] = {0, 0, yShiftIdx, yShiftIdx};
-
-			double locZ[4];
-			locZ[0] = elevation[yIdx*width + xIdx];
-			if (fabs(locZ[0] - no_data) > std::numeric_limits<double>::epsilon())
-			{
-				for (unsigned j=1; j<4; ++j)
-				{
-					locZ[j] = elevation[(yIdx+y_nb[j])*width + (xIdx+x_nb[j])];
-					if (fabs(locZ[j] - no_data) < std::numeric_limits<double>::epsilon())
-						locZ[j]=locZ[0];
-				}
-
-				double ome[4];
-				double xi = 1-fabs(xShift);
-				double eta = 1-fabs(xShift);
-				MathLib::MPhi2D(ome, xi, eta);
-
-				double z(0.0);
-				for(unsigned j = 0; j < 4; ++j)
-					z += ome[j] * locZ[j];
-
-				nodes[i]->updateCoordinates(coords[0], coords[1], z);
-			}
-			else
+		// adjust weights if necessary
+		if (no_data_count > 0)
+		{
+			if (no_data_count == 4) // if there is absolutely no data just use the default value
 			{
-				if (layer_id == 0) // use default value
-					nodes[i]->updateCoordinates(coords[0], coords[1], noDataReplacementValue);
-				else // use z-value from layer above
-					nodes[i]->updateCoordinates(coords[0], coords[1], (*nodes[i-nNodesPerLayer])[2]);
+				nodes[i]->updateCoordinates(coords[0], coords[1], noDataReplacementValue);
+				continue;
 			}
+			const double norm = 4/(4-no_data_count);
+			std::for_each(weight.begin(), weight.end(), [&norm](double &val){val*=norm;});
 		}
 
-		delete raster;
-		return 1;
+		// new value
+		double z = MathLib::scalarProduct<double,4>(weight.data(), pix_val.data());
+		nodes[i]->updateCoordinates(coords[0], coords[1], z);
 	}
-	else
-		ERR("MshLayerMapper::LayerMapping() - Mesh has only %d Layers, cannot assign layer %d.", nLayers, layer_id);
-	return 0;
+
+	delete raster;
+	return 1;
 }
 
-bool MshLayerMapper::isNodeOnRaster(const MeshLib::Node &node,
+bool MeshLayerMapper::isNodeOnRaster(const MeshLib::Node &node,
                                     const std::pair<double, double> &xDim,
                                     const std::pair<double, double> &yDim)
 {
@@ -235,16 +215,16 @@ bool MshLayerMapper::isNodeOnRaster(const MeshLib::Node &node,
 	return true;
 }
 
-MeshLib::Mesh* MshLayerMapper::blendLayersWithSurface(MeshLib::Mesh* mesh, const unsigned nLayers, const std::string &dem_raster)
+MeshLib::Mesh* MeshLayerMapper::blendLayersWithSurface(MeshLib::Mesh &mesh, const unsigned nLayers, const std::string &dem_raster)
 {
 	// construct surface mesh from DEM
 	const MathLib::Vector3 dir(0,0,1);
-	MeshLib::Mesh* dem = MeshLib::MeshSurfaceExtraction::getMeshSurface(*mesh, dir);
-	MshLayerMapper::LayerMapping(dem, dem_raster, 0, 0);
-	const std::vector<MeshLib::Node*> dem_nodes (dem->getNodes());
+	MeshLib::Mesh* dem = MeshLib::MeshSurfaceExtraction::getMeshSurface(mesh, dir);
+	MeshLayerMapper::LayerMapping(*dem, dem_raster, 0, 0);
+	const std::vector<MeshLib::Node*> &dem_nodes (dem->getNodes());
 
-	const std::vector<MeshLib::Node*> mdl_nodes (mesh->getNodes());
-	const size_t nNodes = mesh->getNNodes();
+	const std::vector<MeshLib::Node*> &mdl_nodes (mesh.getNodes());
+	const size_t nNodes (mesh.getNNodes());
 	const size_t nNodesPerLayer = nNodes / (nLayers+1);
 	std::vector<bool> is_surface_node(nNodes, false);
 	std::vector<bool> nodes_below_surface(nNodes, false);
@@ -268,7 +248,7 @@ MeshLib::Mesh* MshLayerMapper::blendLayersWithSurface(MeshLib::Mesh* mesh, const
 	// if node < dem-node: do nothing
 	// if node > dem-node:
 	//		if first node above surface: map to dem and mark as surface node
-	//		else remove node (i.e. don't copy them)
+	//		else remove node (i.e. don't copy it)
 	for (int layer_id=nLayers-1; layer_id>=0; --layer_id)
 	{
 		const size_t firstNode = layer_id * nNodesPerLayer;
@@ -292,6 +272,8 @@ MeshLib::Mesh* MshLayerMapper::blendLayersWithSurface(MeshLib::Mesh* mesh, const
 			}
 		}
 	}
+	delete dem; // no longer needed
+
 	// copy valid nodes to new node vector
 	std::vector<MeshLib::Node*> new_nodes;
 	std::vector<int> node_index_map(nNodes, -1);
@@ -304,8 +286,8 @@ MeshLib::Mesh* MshLayerMapper::blendLayersWithSurface(MeshLib::Mesh* mesh, const
 		}
 
 	// copy elements (old elements need to have at least 4 nodes remaining and form a 3d element
-	const std::vector<MeshLib::Element*> mdl_elements (mesh->getElements());
-	const size_t nElems = mesh->getNElements();
+	const std::vector<MeshLib::Element*> &mdl_elements (mesh.getElements());
+	const size_t nElems (mesh.getNElements());
 	std::vector<MeshLib::Element*> new_elements;
 	for (unsigned j=0; j<nElems; ++j)
 	{
@@ -332,7 +314,7 @@ MeshLib::Mesh* MshLayerMapper::blendLayersWithSurface(MeshLib::Mesh* mesh, const
 				if (!nodes_below_surface[elem->getNode(i)->getID()])
 					top_idx = i-3;
 
-			// construct pyrmid element based on missing node
+			// construct pyramid element based on missing node
 			unsigned idx1 ((top_idx+1)%3);
 			unsigned idx2 ((top_idx+2)%3);
 			e_nodes[0] = new_nodes[node_index_map[elem->getNode(idx1)->getID()]];
diff --git a/Gui/DataView/MshLayerMapper.h b/MeshLib/MeshGenerators/MeshLayerMapper.h
similarity index 81%
rename from Gui/DataView/MshLayerMapper.h
rename to MeshLib/MeshGenerators/MeshLayerMapper.h
index 6f483272ed3e1d47aed1a48bcb098e41ae1ea056..bcbd9a7da94e38a561fa583aec776ed439c7115d 100644
--- a/Gui/DataView/MshLayerMapper.h
+++ b/MeshLib/MeshGenerators/MeshLayerMapper.h
@@ -1,8 +1,8 @@
 /**
- * \file
+ * \file   MeshLayerMapper.h
  * \author Karsten Rink
  * \date   2010-11-01
- * \brief  Definition of the MshLayerMapper class.
+ * \brief  Definition of the MeshLayerMapper class.
  *
  * \copyright
  * Copyright (c) 2013, OpenGeoSys Community (http://www.opengeosys.org)
@@ -12,8 +12,8 @@
  *
  */
 
-#ifndef MSHLAYERMAPPER_H
-#define MSHLAYERMAPPER_H
+#ifndef MESHLAYERMAPPER_H
+#define MESHLAYERMAPPER_H
 
 #include <string>
 
@@ -27,11 +27,11 @@ namespace MeshLib {
 /**
  * \brief Manipulating and adding layers to an existing mesh
  */
-class MshLayerMapper
+class MeshLayerMapper
 {
 public:
-	MshLayerMapper() {}
-	~MshLayerMapper() {}
+	MeshLayerMapper() {}
+	~MeshLayerMapper() {}
 
 	/**
 	 * Based on a triangle-or quad mesh this method creates a 3D mesh with with a given number of prism- or hex-layers
@@ -40,13 +40,13 @@ public:
 	 * \param thickness The thickness of each of these newly added layers
 	 * \return A mesh with the requested number of layers of prism/hex elements
 	 */
-	static MeshLib::Mesh* CreateLayers(const MeshLib::Mesh* mesh, const std::vector<float> &thickness);
+	static MeshLib::Mesh* CreateLayers(const MeshLib::Mesh &mesh, const std::vector<float> &layer_thickness_vector);
 
 	/**
 	 * Maps the z-values of nodes in the designated layer of the given mesh according to the given raster.
 	 * Note: This only results in a valid mesh if the layers don't intersect each other.
 	 */
-	static int LayerMapping(MeshLib::Mesh* msh, const std::string &rasterfile,
+	static int LayerMapping(MeshLib::Mesh &mesh, const std::string &rasterfile,
                             const unsigned nLayers, const unsigned layer_id, double noDataReplacementValue);
 
 	/**
@@ -55,7 +55,7 @@ public:
 	 * remedied at the end of method upon creating the actual mesh from the new node- and element-vector as the mesh-constructor checks for such
 	 * nodes and removes them. This note is just to call this issue to attention in case this methods is changed.
 	 */
-	static MeshLib::Mesh* blendLayersWithSurface(MeshLib::Mesh* mesh, const unsigned nLayers, const std::string &dem_raster);
+	static MeshLib::Mesh* blendLayersWithSurface(MeshLib::Mesh &mesh, const unsigned nLayers, const std::string &dem_raster);
 
 private:
 	/// Checks if the given mesh is within the dimensions given by xDim and yDim.
@@ -64,4 +64,4 @@ private:
 	                           const std::pair<double, double> &yDim);
 };
 
-#endif //MSHLAYERMAPPER_H
+#endif //MESHLAYERMAPPER_H
diff --git a/MeshLib/Node.h b/MeshLib/Node.h
index e6000a6cd5d1e06584a308f48b60168c2a6e40c7..8bdb813c93418f5dc983fde85c54f4038c2f28a4 100644
--- a/MeshLib/Node.h
+++ b/MeshLib/Node.h
@@ -24,9 +24,7 @@
 #include "Mesh.h"
 #include "MeshEditing/removeMeshNodes.h"
 #include "MeshSurfaceExtraction.h"
-#ifdef OGS_BUILD_GUI
-	#include "../Gui/DataView/MshLayerMapper.h"
-#endif
+#include "MeshGenerators/MeshLayerMapper.h"
 
 namespace MeshLib {
 
@@ -38,11 +36,10 @@ class Element;
 class Node : public GeoLib::PointWithID
 {
 	/* friend functions: */
-#ifdef OGS_BUILD_GUI
-	friend int MshLayerMapper::LayerMapping(MeshLib::Mesh* msh, const std::string &rasterfile, const unsigned nLayers,
+	friend int MeshLayerMapper::LayerMapping(MeshLib::Mesh &mesh, const std::string &rasterfile, const unsigned nLayers,
 		                                    const unsigned layer_id, double noDataReplacementValue);
-	friend MeshLib::Mesh* MshLayerMapper::blendLayersWithSurface(MeshLib::Mesh* mesh, const unsigned nLayers, const std::string &dem_raster);
-#endif
+	friend MeshLib::Mesh* MeshLayerMapper::blendLayersWithSurface(MeshLib::Mesh &mesh, const unsigned nLayers, const std::string &dem_raster);
+
 	/* friend classes: */
 	friend class Mesh;
 	friend class MeshRevision;
diff --git a/NumLib/CMakeLists.txt b/NumLib/CMakeLists.txt
index dfa76f13248ea5e4e5ca78d4b98a7bbd1e590768..712d6420aae70374f05ec1b27bf35de3df64fe11 100644
--- a/NumLib/CMakeLists.txt
+++ b/NumLib/CMakeLists.txt
@@ -4,6 +4,14 @@ SET ( SOURCES ${SOURCES_NUMLIB})
 
 GET_SOURCE_FILES(SOURCES_FEM Fem)
 SET ( SOURCES ${SOURCES} ${SOURCES_FEM})
+GET_SOURCE_FILES(SOURCES_FEM_COORDINATESMAPPING Fem/CoordinatesMapping)
+SET ( SOURCES ${SOURCES} ${SOURCES_FEM_COORDINATESMAPPING})
+GET_SOURCE_FILES(SOURCES_FEM_FINTIEELEMENT Fem/FiniteElement)
+SET ( SOURCES ${SOURCES} ${SOURCES_FEM_FINTIEELEMENT})
+GET_SOURCE_FILES(SOURCES_FEM_INTEGRATION Fem/Integration)
+SET ( SOURCES ${SOURCES} ${SOURCES_FEM_INTEGRATION})
+GET_SOURCE_FILES(SOURCES_FEM_SHAPEFUNCTION Fem/ShapeFunction)
+SET ( SOURCES ${SOURCES} ${SOURCES_FEM_SHAPEFUNCTION})
 
 GET_SOURCE_FILES(SOURCES_TIMESTEP TimeStepping)
 GET_SOURCE_FILES(SOURCES_TIMESTEP_ALGORITHMS TimeStepping/Algorithms)
diff --git a/NumLib/Fem/CoordinatesMapping/ShapeMatrices.h b/NumLib/Fem/CoordinatesMapping/ShapeMatrices.h
index c123312cce8dfa57502183946a9a1fc6605e51e4..fb17720d561ce189e8469afa664c2b47ee26d2da 100644
--- a/NumLib/Fem/CoordinatesMapping/ShapeMatrices.h
+++ b/NumLib/Fem/CoordinatesMapping/ShapeMatrices.h
@@ -61,7 +61,10 @@ struct ShapeMatrices
      */
     ShapeMatrices(std::size_t dim, std::size_t n_nodes)
     : N(n_nodes), dNdr(dim, n_nodes), J(dim, dim), detJ(.0),
-      invJ(dim, dim), dNdx(dim, n_nodes) {}
+      invJ(dim, dim), dNdx(dim, n_nodes)
+    {
+        this->setZero();
+    }
 
     ~ShapeMatrices() {}
 
diff --git a/NumLib/Fem/FiniteElement/C0IsoparametricElements.h b/NumLib/Fem/FiniteElement/C0IsoparametricElements.h
index 7635f6ce50aec47d046c9cdeeadf00da17692119..68fb04d38e5a5ea927bc194c9200972f136540be 100644
--- a/NumLib/Fem/FiniteElement/C0IsoparametricElements.h
+++ b/NumLib/Fem/FiniteElement/C0IsoparametricElements.h
@@ -14,8 +14,12 @@
 #ifndef C0ISOPARAMETRICELEMENTS_H_
 #define C0ISOPARAMETRICELEMENTS_H_
 
+#include "MeshLib/Elements/Line.h"
 #include "MeshLib/Elements/Quad.h"
+#include "MeshLib/Elements/Hex.h"
+#include "NumLib/Fem/ShapeFunction/ShapeLine2.h"
 #include "NumLib/Fem/ShapeFunction/ShapeQuad4.h"
+#include "NumLib/Fem/ShapeFunction/ShapeHex8.h"
 #include "NumLib/Fem/Integration/IntegrationGaussRegular.h"
 
 #include "TemplateIsoparametric.h"
@@ -23,12 +27,24 @@
 namespace NumLib
 {
 
+template <class T_SHAPE_VECTOR, class T_DSHAPE_MATRIX, class T_JACOBIAN_MATRIX>
+struct FeLINE2
+{
+    typedef TemplateIsoparametric<MeshLib::Line, ShapeLine2, IntegrationGaussRegular<1>, T_SHAPE_VECTOR, T_DSHAPE_MATRIX, T_JACOBIAN_MATRIX> type;
+};
+
 template <class T_SHAPE_VECTOR, class T_DSHAPE_MATRIX, class T_JACOBIAN_MATRIX>
 struct FeQUAD4
 {
     typedef TemplateIsoparametric<MeshLib::Quad, ShapeQuad4, IntegrationGaussRegular<2>, T_SHAPE_VECTOR, T_DSHAPE_MATRIX, T_JACOBIAN_MATRIX> type;
 };
 
+template <class T_SHAPE_VECTOR, class T_DSHAPE_MATRIX, class T_JACOBIAN_MATRIX>
+struct FeHEX8
+{
+    typedef TemplateIsoparametric<MeshLib::Hex, ShapeHex8, IntegrationGaussRegular<3>, T_SHAPE_VECTOR, T_DSHAPE_MATRIX, T_JACOBIAN_MATRIX> type;
+};
+
 } // NumLib
 
 #endif //C0ISOPARAMETRICELEMENTS_H_
diff --git a/NumLib/Fem/Integration/IntegrationGaussRegular.h b/NumLib/Fem/Integration/IntegrationGaussRegular.h
index 0910cad4156b07b7d3ee0526e741dd2b8b08c421..fd6991ec6368fc83c3619e2d3b6363e53cb0aa32 100644
--- a/NumLib/Fem/Integration/IntegrationGaussRegular.h
+++ b/NumLib/Fem/Integration/IntegrationGaussRegular.h
@@ -49,7 +49,7 @@ public:
     /// Change the integration order.
     void setIntegrationOrder(std::size_t order)
     {
-        this->_n_sampl_pt = std::pow(order, N_DIM);
+        this->_n_sampl_pt = static_cast<std::size_t>(std::pow(order, N_DIM));
         this->_order = order;
     }
 
diff --git a/NumLib/Fem/ShapeFunction/ShapeHex8-impl.h b/NumLib/Fem/ShapeFunction/ShapeHex8-impl.h
new file mode 100644
index 0000000000000000000000000000000000000000..afcf011f759084e1db5da3f343213662b6ea4f9d
--- /dev/null
+++ b/NumLib/Fem/ShapeFunction/ShapeHex8-impl.h
@@ -0,0 +1,61 @@
+/**
+ * \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/project/license
+ *
+ */
+
+namespace NumLib
+{
+
+template <class T_X, class T_N>
+void ShapeHex8::computeShapeFunction(const T_X &r, T_N &N)
+{
+    N[0] = (1.0 - r[0]) * (1.0 - r[1]) * (1.0 - r[2]) * 0.125;
+    N[1] = (1.0 + r[0]) * (1.0 - r[1]) * (1.0 - r[2]) * 0.125;
+    N[2] = (1.0 + r[0]) * (1.0 + r[1]) * (1.0 - r[2]) * 0.125;
+    N[3] = (1.0 - r[0]) * (1.0 + r[1]) * (1.0 - r[2]) * 0.125;
+    N[4] = (1.0 - r[0]) * (1.0 - r[1]) * (1.0 + r[2]) * 0.125;
+    N[5] = (1.0 + r[0]) * (1.0 - r[1]) * (1.0 + r[2]) * 0.125;
+    N[6] = (1.0 + r[0]) * (1.0 + r[1]) * (1.0 + r[2]) * 0.125;
+    N[7] = (1.0 - r[0]) * (1.0 + r[1]) * (1.0 + r[2]) * 0.125;
+}
+
+template <class T_X, class T_N>
+void ShapeHex8::computeGradShapeFunction(const T_X &r, T_N &dN)
+{
+    // dN/dx
+    dN[0] = -(1.0 - r[1]) * (1.0 - r[2]) * 0.125;
+    dN[1] = -dN[0];
+    dN[2] = +(1.0 + r[1]) * (1.0 - r[2]) * 0.125;
+    dN[3] = -dN[2];
+    dN[4] = -(1.0 - r[1]) * (1.0 + r[2]) * 0.125;
+    dN[5] = -dN[4];
+    dN[6] = +(1.0 + r[1]) * (1.0 + r[2]) * 0.125;
+    dN[7] = -dN[6];
+
+    // dN/dy
+    dN[8]  = -(1.0 - r[0]) * (1.0 - r[2]) * 0.125;
+    dN[9]  = -(1.0 + r[0]) * (1.0 - r[2]) * 0.125;
+    dN[10] = -dN[9];
+    dN[11] = -dN[8];
+    dN[12] = -(1.0 - r[0]) * (1.0 + r[2]) * 0.125;
+    dN[13] = -(1.0 + r[0]) * (1.0 + r[2]) * 0.125;
+    dN[14] = -dN[13];
+    dN[15] = -dN[12];
+
+    // dN/dz
+    dN[16] = -(1.0 - r[0]) * (1.0 - r[1]) * 0.125;
+    dN[17] = -(1.0 + r[0]) * (1.0 - r[1]) * 0.125;
+    dN[18] = -(1.0 + r[0]) * (1.0 + r[1]) * 0.125;
+    dN[19] = -(1.0 - r[0]) * (1.0 + r[1]) * 0.125;
+    dN[20] = -dN[16];
+    dN[21] = -dN[17];
+    dN[22] = -dN[18];
+    dN[23] = -dN[19];
+}
+
+}
+
diff --git a/NumLib/Fem/ShapeFunction/ShapeHex8.h b/NumLib/Fem/ShapeFunction/ShapeHex8.h
new file mode 100644
index 0000000000000000000000000000000000000000..5a5d12bd9e91f0c09fd65d53e1556e40cb38c29c
--- /dev/null
+++ b/NumLib/Fem/ShapeFunction/ShapeHex8.h
@@ -0,0 +1,64 @@
+/**
+ * \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/project/license
+ *
+ */
+
+#ifndef SHAPEHEX8_H_
+#define SHAPEHEX8_H_
+
+namespace NumLib
+{
+
+/**
+ *  Shape function for a 8-nodes hex element in natural coordinates
+ *
+ * \verbatim
+ *       (-1, 1, 1) 7-----------6 ( 1, 1, 1)
+ *                 /:          /|
+ *                / :         / |
+ *               /  :        /  |
+ *              /   :       /   |
+ *             /    :      /    |
+ * (-1,-1, 1) 4-----------5 ( 1,-1, 1)
+ *            |     :     |     |
+ *       (-1, 1,-1) 3.....|.....2 ( 1, 1,-1)
+ *            |    .      |    /
+ *            |   .       |   /
+ *            |  .        |  /
+ *            | .         | /
+ *            |.          |/
+ * (-1,-1,-1) 0-----------1 ( 1,-1,-1)
+ *          0
+ * \endverbatim
+ */
+class ShapeHex8
+{
+public:
+    /**
+     * Evaluate the shape function at the given point
+     *
+     * @param [in]  r   natural coordinates (r,s,t)
+     * @param [out] N   a vector of calculated shape functions
+     */
+    template <class T_X, class T_N>
+    static void computeShapeFunction(const T_X &r, T_N &N);
+
+    /**
+     * Evaluate derivatives of the shape function at the given point
+     *
+     * @param [in]  r   natural coordinates (r,s,t)
+     * @param [out] dN  a matrix of the derivatives
+     */
+    template <class T_X, class T_N>
+    static void computeGradShapeFunction(const T_X &r, T_N &dN);
+};
+
+}
+
+#include "ShapeHex8-impl.h"
+
+#endif //SHAPEHEX8_H_
diff --git a/NumLib/Fem/ShapeFunction/ShapeLine2-impl.h b/NumLib/Fem/ShapeFunction/ShapeLine2-impl.h
new file mode 100644
index 0000000000000000000000000000000000000000..cc458489f10792d30cb7f14a02d8eabe2e35248a
--- /dev/null
+++ b/NumLib/Fem/ShapeFunction/ShapeLine2-impl.h
@@ -0,0 +1,28 @@
+/**
+ * \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/project/license
+ *
+ */
+
+namespace NumLib
+{
+
+template <class T_X, class T_N>
+void ShapeLine2::computeShapeFunction(const T_X &r, T_N &N)
+{
+    N[0] = (1.0 - r[0]) * 0.5;
+    N[1] = (1.0 + r[0]) * 0.5;
+}
+
+template <class T_X, class T_N>
+void ShapeLine2::computeGradShapeFunction(const T_X &/*r*/, T_N &dN)
+{
+    dN[0] = -0.5;
+    dN[1] = 0.5;
+}
+
+}
+
diff --git a/NumLib/Fem/ShapeFunction/ShapeLine2.h b/NumLib/Fem/ShapeFunction/ShapeLine2.h
new file mode 100644
index 0000000000000000000000000000000000000000..13b70e6c299f75cae058bac3cd4ec7a9da953094
--- /dev/null
+++ b/NumLib/Fem/ShapeFunction/ShapeLine2.h
@@ -0,0 +1,50 @@
+/**
+ * \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/project/license
+ *
+ */
+
+#ifndef SHAPELINE2_H_
+#define SHAPELINE2_H_
+
+namespace NumLib
+{
+
+/**
+ *  Shape function for a line element of two nodes in natural coordinates
+ *
+ * \verbatim
+ *  2 (-1,0)     1 (1,0)
+ *     *--------*
+ * \endverbatim
+ */
+class ShapeLine2
+{
+public:
+    /**
+     * Evaluate the shape function at the given point
+     *
+     * @param [in]  r    point coordinates
+     * @param [out] N   a vector of calculated shape function.
+     */
+    template <class T_X, class T_N>
+    static void computeShapeFunction(const T_X &r, T_N &N);
+
+    /**
+     * Evaluate derivatives of the shape function at the given point
+     *
+     * @param [in]  r    point coordinates
+     * @param [out] dN  a matrix of the derivatives
+     */
+    template <class T_X, class T_N>
+    static void computeGradShapeFunction(const T_X &r, T_N &dN);
+};
+
+}
+
+#include "ShapeLine2-impl.h"
+
+#endif //SHAPELINE2_H_
diff --git a/Tests/MeshLib/TestTriLineMesh.cpp b/Tests/MeshLib/TestTriLineMesh.cpp
index ce9997241e6e9a2ec01ecb6bcb76eb17214ca04b..a745bc8779f373015c33deff62f69134ce63d670 100644
--- a/Tests/MeshLib/TestTriLineMesh.cpp
+++ b/Tests/MeshLib/TestTriLineMesh.cpp
@@ -73,12 +73,6 @@ class MeshLibTriLineMesh : public ::testing::Test
             elements[e])
                 != nodes[n]->getElements().cend();
     }
-
-    bool
-    elementIsNeigborOf(std::size_t const a, std::size_t const b) const
-    {
-        return false;
-    }
 };
 
 TEST_F(MeshLibTriLineMesh, Construction)
diff --git a/Tests/NumLib/CoordinatesMappingTestData/TestHex8.h b/Tests/NumLib/CoordinatesMappingTestData/TestHex8.h
new file mode 100644
index 0000000000000000000000000000000000000000..c164026df6e2508b622c99d975d9925e241f9cd9
--- /dev/null
+++ b/Tests/NumLib/CoordinatesMappingTestData/TestHex8.h
@@ -0,0 +1,126 @@
+/**
+ * \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
+ */
+
+#ifndef COORDINATESMAPPINGTESTDATA_TESTHEX8_H_
+#define COORDINATESMAPPINGTESTDATA_TESTHEX8_H_
+
+#include "MeshLib/Elements/Hex.h"
+#include "NumLib/Fem/ShapeFunction/ShapeHex8.h"
+
+namespace CoordinatesMappingTestData
+{
+
+class TestHex8
+{
+ public:
+    // Element information
+    typedef MeshLib::Hex ElementType;
+    typedef NumLib::ShapeHex8 ShapeFunctionType;
+    static const unsigned dim = 3; //ElementType::dimension;
+    static const unsigned e_nnodes = ElementType::n_all_nodes;
+    // Coordinates where shape functions are evaluated
+    static const double r[dim];
+    // Expected results for natural shape
+    static const double nat_exp_N[e_nnodes];
+    static const double nat_exp_dNdr[e_nnodes*dim];
+    // Expected results for irregular shape
+    static const double ir_exp_J[dim*dim];
+    static const double ir_exp_invJ[dim*dim];
+    static const double ir_exp_detJ;
+    static const double ir_exp_dNdx[e_nnodes*dim];
+    // Expected results for clock-wise node ordering
+    static const double cl_exp_J[dim*dim];
+    // Expected results for zero volume
+    static const double cl_exp_detJ;
+    static const double ze_exp_J[dim*dim];
+
+    // element shape identical to that in natural coordinates (see ShapeHex8.h)
+    MeshLib::Hex* createNaturalShape()
+    {
+        MeshLib::Node** nodes = new MeshLib::Node*[e_nnodes];
+        nodes[0] = new MeshLib::Node(-1.0, -1.0, -1.0);
+        nodes[1] = new MeshLib::Node( 1.0, -1.0, -1.0);
+        nodes[2] = new MeshLib::Node( 1.0,  1.0, -1.0);
+        nodes[3] = new MeshLib::Node(-1.0,  1.0, -1.0);
+        nodes[4] = new MeshLib::Node(-1.0, -1.0,  1.0);
+        nodes[5] = new MeshLib::Node( 1.0, -1.0,  1.0);
+        nodes[6] = new MeshLib::Node( 1.0,  1.0,  1.0);
+        nodes[7] = new MeshLib::Node(-1.0,  1.0,  1.0);
+        return new MeshLib::Hex(nodes);
+    }
+
+    // element having irregular or skew shape
+    MeshLib::Hex* createIrregularShape()
+    {
+        // two times longer in z direction than the natural
+        MeshLib::Node** nodes = new MeshLib::Node*[e_nnodes];
+        nodes[0] = new MeshLib::Node(-1.0, -1.0, -1.0);
+        nodes[1] = new MeshLib::Node( 1.0, -1.0, -1.0);
+        nodes[2] = new MeshLib::Node( 1.0,  1.0, -1.0);
+        nodes[3] = new MeshLib::Node(-1.0,  1.0, -1.0);
+        nodes[4] = new MeshLib::Node(-1.0, -1.0,  3.0);
+        nodes[5] = new MeshLib::Node( 1.0, -1.0,  3.0);
+        nodes[6] = new MeshLib::Node( 1.0,  1.0,  3.0);
+        nodes[7] = new MeshLib::Node(-1.0,  1.0,  3.0);
+        return new MeshLib::Hex(nodes);
+    }
+
+    // invalid case: clock wise node ordering
+    MeshLib::Hex* createClockWise()
+    {
+        MeshLib::Node** nodes = new MeshLib::Node*[e_nnodes];
+        nodes[0] = new MeshLib::Node(-1.0, -1.0, -1.0);
+        nodes[3] = new MeshLib::Node( 1.0, -1.0, -1.0);
+        nodes[2] = new MeshLib::Node( 1.0,  1.0, -1.0);
+        nodes[1] = new MeshLib::Node(-1.0,  1.0, -1.0);
+        nodes[4] = new MeshLib::Node(-1.0, -1.0,  1.0);
+        nodes[7] = new MeshLib::Node( 1.0, -1.0,  1.0);
+        nodes[6] = new MeshLib::Node( 1.0,  1.0,  1.0);
+        nodes[5] = new MeshLib::Node(-1.0,  1.0,  1.0);
+        return new MeshLib::Hex(nodes);
+    }
+
+    // invalid case: zero volume
+    MeshLib::Hex* createZeroVolume()
+    {
+        MeshLib::Node** nodes = new MeshLib::Node*[e_nnodes];
+        nodes[0] = new MeshLib::Node(-1.0, -1.0,  1.0);
+        nodes[1] = new MeshLib::Node( 1.0, -1.0,  1.0);
+        nodes[2] = new MeshLib::Node( 1.0,  1.0,  1.0);
+        nodes[3] = new MeshLib::Node(-1.0,  1.0,  1.0);
+        nodes[4] = new MeshLib::Node(-1.0, -1.0,  1.0);
+        nodes[5] = new MeshLib::Node( 1.0, -1.0,  1.0);
+        nodes[6] = new MeshLib::Node( 1.0,  1.0,  1.0);
+        nodes[7] = new MeshLib::Node(-1.0,  1.0,  1.0);
+        return new MeshLib::Hex(nodes);
+    }
+};
+
+const double TestHex8::r[dim] = {0.5, 0.5, 0.5};
+const double TestHex8::nat_exp_N[e_nnodes]
+    = {0.015625, 0.046875, 0.140625, 0.046875, 0.046875, 0.140625, 0.421875, 0.140625};
+const double TestHex8::nat_exp_dNdr[e_nnodes*dim]
+    = {-0.03125, 0.03125, 0.09375, -0.09375, -0.09375, 0.09375, 0.28125, -0.28125,
+       -0.03125, -0.09375, 0.09375, 0.03125, -0.09375, -0.28125,  0.28125,  0.09375,
+       -0.03125, -0.09375, -0.28125, -0.09375,  0.03125,  0.09375,  0.28125,  0.09375};
+const double TestHex8::ir_exp_J[dim*dim]= {1.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 2.0};
+const double TestHex8::ir_exp_detJ = 2.;
+const double TestHex8::ir_exp_invJ[dim*dim] = {1.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 1./2.0};
+const double TestHex8::ir_exp_dNdx[dim*e_nnodes]
+    = {-0.03125, 0.03125, 0.09375, -0.09375, -0.09375, 0.09375, 0.28125, -0.28125,
+       -0.03125, -0.09375, 0.09375, 0.03125, -0.09375, -0.28125,  0.28125,  0.09375,
+       -0.015625, -0.046875, -0.140625, -0.046875,  0.015625,  0.046875,  0.140625,  0.046875};
+const double TestHex8::cl_exp_J[dim*dim] = {0.0, 1.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 1.0};
+const double TestHex8::cl_exp_detJ = -1.;
+const double TestHex8::ze_exp_J[dim*dim] = {1.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0};
+
+}
+
+#endif //OGS_USE_EIGEN
+
+
diff --git a/Tests/NumLib/CoordinatesMappingTestData/TestLine2.h b/Tests/NumLib/CoordinatesMappingTestData/TestLine2.h
new file mode 100644
index 0000000000000000000000000000000000000000..048380919c5ca7d314a9c8088aec266ab6a44829
--- /dev/null
+++ b/Tests/NumLib/CoordinatesMappingTestData/TestLine2.h
@@ -0,0 +1,93 @@
+/**
+ * \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
+ */
+
+#ifndef COORDINATESMAPPINGTESTDATA_TESTLINE2_H_
+#define COORDINATESMAPPINGTESTDATA_TESTLINE2_H_
+
+#include "MeshLib/Elements/Line.h"
+#include "NumLib/Fem/ShapeFunction/ShapeLine2.h"
+
+namespace CoordinatesMappingTestData
+{
+
+class TestLine2
+{
+public:
+    // Element information
+    typedef MeshLib::Line ElementType;
+    typedef NumLib::ShapeLine2 ShapeFunctionType;
+    static const unsigned dim = ElementType::dimension;
+    static const unsigned e_nnodes = ElementType::n_all_nodes;
+    // Coordinates where shape functions are evaluated
+    static const double r[dim];
+    // Expected results for natural shape
+    static const double nat_exp_N[e_nnodes];
+    static const double nat_exp_dNdr[e_nnodes*dim];
+    // Expected results for irregular shape
+    static const double ir_exp_J[dim*dim];
+    static const double ir_exp_invJ[dim*dim];
+    static const double ir_exp_detJ;
+    static const double ir_exp_dNdx[e_nnodes*dim];
+    // Expected results for clock-wise node ordering
+    static const double cl_exp_J[dim*dim];
+    // Expected results for zero volume
+    static const double cl_exp_detJ;
+    static const double ze_exp_J[dim*dim];
+
+    // element shape identical to that in natural coordinates (see ShapeLine2.h)
+    static MeshLib::Line* createNaturalShape()
+    {
+        MeshLib::Node** nodes = new MeshLib::Node*[e_nnodes];
+        nodes[0] = new MeshLib::Node(-1.0, 0.0, 0.0);
+        nodes[1] = new MeshLib::Node( 1.0, 0.0, 0.0);
+        return new MeshLib::Line(nodes);
+    }
+
+    // element having irregular or skew shape
+    MeshLib::Line* createIrregularShape()
+    {
+        // two times longer than the natural
+        MeshLib::Node** nodes = new MeshLib::Node*[e_nnodes];
+        nodes[0] = new MeshLib::Node(-2.0, 0.0, 0.0);
+        nodes[1] = new MeshLib::Node( 2.0, 0.0, 0.0);
+        return new MeshLib::Line(nodes);
+    }
+
+    // invalid case: clock wise node ordering
+    MeshLib::Line* createClockWise()
+    {
+        MeshLib::Node** nodes = new MeshLib::Node*[e_nnodes];
+        nodes[1] = new MeshLib::Node(-1.0, 0.0, 0.0);
+        nodes[0] = new MeshLib::Node( 1.0, 0.0, 0.0);
+        return new MeshLib::Line(nodes);
+    }
+
+    // invalid case: zero volume
+    MeshLib::Line* createZeroVolume()
+    {
+        MeshLib::Node** nodes = new MeshLib::Node*[e_nnodes];
+        nodes[0] = new MeshLib::Node(0.0, 0.0, 0.0);
+        nodes[1] = new MeshLib::Node(0.0, 0.0, 0.0);
+        return new MeshLib::Line(nodes);
+    }
+};
+
+const double TestLine2::r[dim] = {0.5};
+const double TestLine2::nat_exp_N[e_nnodes] = {0.25, 0.75};
+const double TestLine2::nat_exp_dNdr[e_nnodes*dim] = {-0.5, 0.5};
+const double TestLine2::ir_exp_J[dim*dim] = {2.0};
+const double TestLine2::ir_exp_invJ[dim*dim] = {0.5};
+const double TestLine2::ir_exp_detJ = 2.;
+const double TestLine2::ir_exp_dNdx[dim*e_nnodes] = {-0.25, 0.25};
+const double TestLine2::cl_exp_J[dim*dim] = {-1.};
+const double TestLine2::cl_exp_detJ = -1;
+const double TestLine2::ze_exp_J[dim*dim] = {0.0};
+
+}
+
+#endif
diff --git a/Tests/NumLib/CoordinatesMappingTestData/TestQuad4.h b/Tests/NumLib/CoordinatesMappingTestData/TestQuad4.h
new file mode 100644
index 0000000000000000000000000000000000000000..b2548cab02fccca6ea911c86ca8327cd565c8327
--- /dev/null
+++ b/Tests/NumLib/CoordinatesMappingTestData/TestQuad4.h
@@ -0,0 +1,105 @@
+/**
+ * \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
+ */
+
+#ifndef COORDINATESMAPPINGTESTDATA_TESTQUAD4_H_
+#define COORDINATESMAPPINGTESTDATA_TESTQUAD4_H_
+
+#include "MeshLib/Elements/Quad.h"
+#include "NumLib/Fem/ShapeFunction/ShapeQuad4.h"
+
+namespace CoordinatesMappingTestData
+{
+
+class TestQuad4
+{
+ public:
+    // Element information
+    typedef MeshLib::Quad ElementType;
+    typedef NumLib::ShapeQuad4 ShapeFunctionType;
+    static const unsigned dim = 2; //ElementType::dimension;
+    static const unsigned e_nnodes = ElementType::n_all_nodes;
+    // Coordinates where shape functions are evaluated
+    static const double r[dim];
+    // Expected results for natural shape
+    static const double nat_exp_N[e_nnodes];
+    static const double nat_exp_dNdr[e_nnodes*dim];
+    // Expected results for irregular shape
+    static const double ir_exp_J[dim*dim];
+    static const double ir_exp_invJ[dim*dim];
+    static const double ir_exp_detJ;
+    static const double ir_exp_dNdx[e_nnodes*dim];
+    // Expected results for clock-wise node ordering
+    static const double cl_exp_J[dim*dim];
+    // Expected results for zero volume
+    static const double cl_exp_detJ;
+    static const double ze_exp_J[dim*dim];
+
+    // element shape identical to that in natural coordinates
+    MeshLib::Quad* createNaturalShape()
+    {
+        MeshLib::Node** nodes = new MeshLib::Node*[e_nnodes];
+        nodes[0] = new MeshLib::Node( 1.0,  1.0,  0.0);
+        nodes[1] = new MeshLib::Node(-1.0,  1.0,  0.0);
+        nodes[2] = new MeshLib::Node(-1.0, -1.0,  0.0);
+        nodes[3] = new MeshLib::Node( 1.0, -1.0,  0.0);
+        return new MeshLib::Quad(nodes);
+    }
+
+    // element having irregular or skew shape
+    MeshLib::Quad* createIrregularShape()
+    {
+        MeshLib::Node** nodes = new MeshLib::Node*[e_nnodes];
+        nodes[0] = new MeshLib::Node(-0.5, -0.5,  0.0);
+        nodes[1] = new MeshLib::Node( 0.6, -0.6,  0.0);
+        nodes[2] = new MeshLib::Node( 0.5,  0.4,  0.0);
+        nodes[3] = new MeshLib::Node(-0.3,  0.1,  0.0);
+        return new MeshLib::Quad(nodes);
+    }
+
+    // invalid case: clock wise node ordering
+    MeshLib::Quad* createClockWise()
+    {
+        MeshLib::Node** nodes = new MeshLib::Node*[e_nnodes];
+        nodes[0] = new MeshLib::Node( 1.0,  1.0,  0.0);
+        nodes[3] = new MeshLib::Node(-1.0,  1.0,  0.0);
+        nodes[2] = new MeshLib::Node(-1.0, -1.0,  0.0);
+        nodes[1] = new MeshLib::Node( 1.0, -1.0,  0.0);
+        return new MeshLib::Quad(nodes);
+    }
+
+    // invalid case: zero area
+    MeshLib::Quad* createZeroVolume()
+    {
+        MeshLib::Node** nodes = new MeshLib::Node*[e_nnodes];
+        nodes[0] = new MeshLib::Node( 1.0,  1.0,  0.0);
+        nodes[1] = new MeshLib::Node(-1.0,  1.0,  0.0);
+        nodes[2] = new MeshLib::Node(-1.0,  1.0,  0.0);
+        nodes[3] = new MeshLib::Node( 1.0,  1.0,  0.0);
+        return new MeshLib::Quad(nodes);
+    }
+};
+
+const double TestQuad4::r[dim] = {0.5, 0.5};
+const double TestQuad4::nat_exp_N[e_nnodes] = {0.5625, 0.1875, 0.0625, 0.1875};
+const double TestQuad4::nat_exp_dNdr[e_nnodes*dim]
+    = {0.375, -0.375, -0.125, 0.125, 0.375, 0.125, -0.125, -0.375};
+const double TestQuad4::ir_exp_J[dim*dim] = {-0.5125, 0.0, -0.0625, -0.35};
+const double TestQuad4::ir_exp_detJ = 0.179375;
+const double TestQuad4::ir_exp_invJ[dim*dim]
+    = {-1.9512195121951219, 0.0, 0.3484320557491290, -2.8571428571428572};
+const double TestQuad4::ir_exp_dNdx[dim*e_nnodes]
+    = {-0.73170731707317072, 0.73170731707317072, 0.243902439024390, -0.24390243902439029,
+       -0.940766550522648, -0.48780487804878048, 0.313588850174216, 1.1149825783972125};
+const double TestQuad4::cl_exp_J[dim*dim] = {0.0, 1.0, 1.0, 0.0};
+const double TestQuad4::cl_exp_detJ = -1.;
+const double TestQuad4::ze_exp_J[dim*dim] = {1.0, 0.0, 0.0, 0.0};
+
+}
+
+#endif
+
diff --git a/Tests/NumLib/FeTestData/MatrixTools.h b/Tests/NumLib/FeTestData/MatrixTools.h
new file mode 100644
index 0000000000000000000000000000000000000000..b4f986c09c6856b1d68bf43303708205d6cef5ff
--- /dev/null
+++ b/Tests/NumLib/FeTestData/MatrixTools.h
@@ -0,0 +1,39 @@
+/**
+ * \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/project/license
+ *
+ */
+
+#ifndef MATRIXTOOLS_H_
+#define MATRIXTOOLS_H_
+
+namespace FeTestData
+{
+
+// copy matrix entries in upper triangle to lower triangle
+template <class T_MATRIX, typename ID_TYPE=signed>
+inline void copyUpperToLower(const ID_TYPE dim, T_MATRIX &m)
+{
+    for (ID_TYPE i=0; i<dim; i++)
+        for (ID_TYPE j=0; j<i; j++)
+            m(i,j) = m(j,i);
+}
+
+// set an identity matrix
+template <class T_MATRIX, typename ID_TYPE=signed>
+inline void setIdentityMatrix(unsigned dim, T_MATRIX &m)
+{
+    for (unsigned i=0; i<dim; i++)
+        for (unsigned j=0; j<dim; j++)
+            m(i,j) = 0.0;
+    for (unsigned i=0; i<dim; i++)
+        m(i,i) = 1.0;
+}
+
+} // namespace
+
+#endif
+
diff --git a/Tests/NumLib/FeTestData/TestFeHEX8.h b/Tests/NumLib/FeTestData/TestFeHEX8.h
new file mode 100644
index 0000000000000000000000000000000000000000..e566087f9864633e5f98603f09d544a0d6c2b8fe
--- /dev/null
+++ b/Tests/NumLib/FeTestData/TestFeHEX8.h
@@ -0,0 +1,93 @@
+/**
+ * \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/project/license
+ *
+ */
+
+#ifndef TESTFEHEX8_H_
+#define TESTFEHEX8_H_
+
+#include "MeshLib/Elements/Hex.h"
+#include "NumLib/Fem/FiniteElement/C0IsoparametricElements.h"
+
+#include "MatrixTools.h"
+
+namespace FeTestData
+{
+
+// Test case for HEX8
+class TestFeHEX8
+{
+public:
+    // Fe type information
+    template <class T_MATRIX_TYPES>
+    struct FeType
+    {
+        typedef NumLib::FeHEX8<typename T_MATRIX_TYPES::NodalVectorType, typename T_MATRIX_TYPES::DimNodalMatrixType, typename T_MATRIX_TYPES::DimMatrixType> type;
+    };
+    typedef MeshLib::Hex MeshElementType;
+    static const unsigned dim = 3; //MeshElementType::dimension;
+    static const unsigned e_nnodes = MeshElementType::n_all_nodes;
+    static const unsigned n_sample_pt_order2 = 2*2*2;
+    static const unsigned n_sample_pt_order3 = 3*3*3;
+
+    /// create a mesh element
+    MeshLib::Hex* createMeshElement()
+    {
+        // cubic
+        MeshLib::Node** nodes = new MeshLib::Node*[e_nnodes];
+        nodes[0] = new MeshLib::Node(0.0, 0.0, 0.0);
+        nodes[1] = new MeshLib::Node(1.0, 0.0, 0.0);
+        nodes[2] = new MeshLib::Node(1.0, 1.0, 0.0);
+        nodes[3] = new MeshLib::Node(0.0, 1.0, 0.0);
+        nodes[4] = new MeshLib::Node(0.0, 0.0, 1.0);
+        nodes[5] = new MeshLib::Node(1.0, 0.0, 1.0);
+        nodes[6] = new MeshLib::Node(1.0, 1.0, 1.0);
+        nodes[7] = new MeshLib::Node(0.0, 1.0, 1.0);
+        return new MeshLib::Hex(nodes);
+    }
+
+    /// set an expected mass matrix
+    template <class T_MATRIX>
+    void setExpectedMassMatrix(T_MATRIX &m)
+    {
+        // set upper triangle entries
+        m(0,0) = 1.0; m(0,1) = 1./2; m(0,2) = 1./4; m(0,3) = 1./2; m(0,4) = 1./2; m(0,5) = 1./4; m(0,6) = 1./8; m(0,7) = 1./4;
+        m(1,1) = 1.0; m(1,2) = 1./2; m(1,3) = 1./4; m(1,4) = 1./4; m(1,5) = 1./2; m(1,6) = 1./4; m(1,7) = 1./8;
+        m(2,2) = 1.0; m(2,3) = 1./2; m(2,4) = 1./8; m(2,5) = 1./4; m(2,6) = 1./2; m(2,7) = 1./4;
+        m(3,3) = 1.0; m(3,4) = 1./4; m(3,5) = 1./8; m(3,6) = 1./4; m(3,7) = 1./2;
+        m(4,4) = 1.0; m(4,5) = 1./2; m(4,6) = 1./4; m(4,7) = 1./2;
+        m(5,5) = 1.0; m(5,6) = 1./2; m(5,7) = 1./4;
+        m(6,6) = 1.0; m(6,7) = 1./2;
+        m(7,7) = 1.0;
+        // make symmetric
+        copyUpperToLower(e_nnodes, m);
+        m *= 1./27.;
+    }
+
+    /// set an expected laplace matrix
+    template <class T_MATRIX>
+    void setExpectedLaplaceMatrix(double k, T_MATRIX &m)
+    {
+        // set upper triangle entries
+        m(0,0) = 2./3; m(0,1) = 0; m(0,2) = -1./6; m(0,3) = 0; m(0,4) = 0; m(0,5) = -1./6; m(0,6) = -1./6; m(0,7) = -1./6;
+        m(1,1) = 2./3; m(1,2) = 0; m(1,3) = -1./6; m(1,4) = -1./6; m(1,5) = 0; m(1,6) = -1./6; m(1,7) = -1./6;
+        m(2,2) = 2./3; m(2,3) = 0; m(2,4) = -1./6; m(2,5) = -1./6; m(2,6) = 0; m(2,7) = -1./6;
+        m(3,3) = 2./3; m(3,4) = -1./6; m(3,5) = -1./6; m(3,6) = -1./6; m(3,7) = 0;
+        m(4,4) = 2./3; m(4,5) = 0; m(4,6) = -1./6; m(4,7) = 0;
+        m(5,5) = 2./3; m(5,6) = 0; m(5,7) = -1./6;
+        m(6,6) = 2./3; m(6,7) = 0;
+        m(7,7) = 2./3;
+        // make symmetric
+        copyUpperToLower(e_nnodes, m);
+        m *= k/2.;
+    }
+};
+
+} // namespace
+
+#endif
+
diff --git a/Tests/NumLib/FeTestData/TestFeLINE2.h b/Tests/NumLib/FeTestData/TestFeLINE2.h
new file mode 100644
index 0000000000000000000000000000000000000000..ce7b5e68d36dc79822a1591ace1c8211844cc982
--- /dev/null
+++ b/Tests/NumLib/FeTestData/TestFeLINE2.h
@@ -0,0 +1,73 @@
+/**
+ * \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/project/license
+ *
+ */
+
+#ifndef TESTFELINE2_H_
+#define TESTFELINE2_H_
+
+#include "MeshLib/Elements/Line.h"
+#include "NumLib/Fem/FiniteElement/C0IsoparametricElements.h"
+
+#include "MatrixTools.h"
+
+namespace FeTestData
+{
+
+// Test case for LINE2
+class TestFeLINE2
+{
+public:
+    // Fe type information
+    template <class T_MATRIX_TYPES>
+    struct FeType
+    {
+        typedef NumLib::FeLINE2<typename T_MATRIX_TYPES::NodalVectorType, typename T_MATRIX_TYPES::DimNodalMatrixType, typename T_MATRIX_TYPES::DimMatrixType> type;
+    };
+    typedef MeshLib::Line MeshElementType;
+    static const unsigned dim = MeshElementType::dimension;
+    static const unsigned e_nnodes = MeshElementType::n_all_nodes;
+    static const unsigned n_sample_pt_order2 = 2;
+    static const unsigned n_sample_pt_order3 = 3;
+
+    /// create a mesh element
+    MeshLib::Line* createMeshElement()
+    {
+        MeshLib::Node** nodes = new MeshLib::Node*[e_nnodes];
+        nodes[0] = new MeshLib::Node(0.0, 0.0, 0.0);
+        nodes[1] = new MeshLib::Node(1.0, 0.0, 0.0);
+        return new MeshLib::Line(nodes);
+    }
+
+    /// set an expected mass matrix
+    template <class T_MATRIX>
+    void setExpectedMassMatrix(T_MATRIX &m)
+    {
+        // set upper triangle entries
+        m(0,0) = 1./3.; m(0,1) = 1./6.;
+        m(1,1) = 1./3.;
+        // make symmetric
+        copyUpperToLower(2, m);
+    }
+
+    /// set an expected laplace matrix
+    template <class T_MATRIX>
+    void setExpectedLaplaceMatrix(double k, T_MATRIX &m)
+    {
+        // set upper triangle entries
+        m(0,0) = 1.0; m(0,1) = -1.0;
+        m(1,1) = 1.0;
+        // make symmetric
+        copyUpperToLower(2, m);
+        m *= k;
+    }
+};
+
+} // namespace
+
+#endif
+
diff --git a/Tests/NumLib/FeTestData/TestFeQUAD4.h b/Tests/NumLib/FeTestData/TestFeQUAD4.h
new file mode 100644
index 0000000000000000000000000000000000000000..009e531ad0afd3d1908240a683ec05410da49598
--- /dev/null
+++ b/Tests/NumLib/FeTestData/TestFeQUAD4.h
@@ -0,0 +1,81 @@
+/**
+ * \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/project/license
+ *
+ */
+
+#ifndef TESTFEQUAD4_H_
+#define TESTFEQUAD4_H_
+
+#include "MeshLib/Elements/Quad.h"
+#include "NumLib/Fem/FiniteElement/C0IsoparametricElements.h"
+
+#include "MatrixTools.h"
+
+namespace FeTestData
+{
+
+// Test case for QUAD4
+class TestFeQUAD4
+{
+public:
+    // Fe type information
+    template <class T_MATRIX_TYPES>
+    struct FeType
+    {
+        typedef NumLib::FeQUAD4<typename T_MATRIX_TYPES::NodalVectorType, typename T_MATRIX_TYPES::DimNodalMatrixType, typename T_MATRIX_TYPES::DimMatrixType> type;
+    };
+    typedef MeshLib::Quad MeshElementType;
+    static const unsigned dim = 2; //MeshElementType::dimension;
+    static const unsigned e_nnodes = MeshElementType::n_all_nodes;
+    static const unsigned n_sample_pt_order2 = 2*2;
+    static const unsigned n_sample_pt_order3 = 3*3;
+
+    /// create a mesh element
+    MeshLib::Quad* createMeshElement()
+    {
+        // square
+        MeshLib::Node** nodes = new MeshLib::Node*[e_nnodes];
+        nodes[0] = new MeshLib::Node(0.0, 0.0, 0.0);
+        nodes[1] = new MeshLib::Node(1.0, 0.0, 0.0);
+        nodes[2] = new MeshLib::Node(1.0, 1.0, 0.0);
+        nodes[3] = new MeshLib::Node(0.0, 1.0, 0.0);
+        return new MeshLib::Quad(nodes);
+    }
+
+    /// set an expected mass matrix for 1m x 1m
+    template <class T_MATRIX>
+    void setExpectedMassMatrix(T_MATRIX &m)
+    {
+        // set upper triangle entries
+        m(0,0) = 1.0; m(0,1) = 1./2; m(0,2) = 1./4; m(0,3) = 1./2;
+        m(1,1) = 1.0; m(1,2) = 1./2; m(1,3) = 1./4;
+        m(2,2) = 1.0; m(2,3) = 1./2;
+        m(3,3) = 1.0;
+        // make symmetric
+        copyUpperToLower(4, m);
+        m *= 1./9.;
+    }
+
+    /// set an expected laplace matrix for 1m x 1m
+    template <class T_MATRIX>
+    void setExpectedLaplaceMatrix(double k, T_MATRIX &m)
+    {
+        // set upper triangle entries
+        m(0,0) = 4.0; m(0,1) = -1.0; m(0,2) = -2.0; m(0,3) = -1.0;
+        m(1,1) = 4.0; m(1,2) = -1.0; m(1,3) = -2.0;
+        m(2,2) = 4.0; m(2,3) = -1.0;
+        m(3,3) = 4.0;
+        // make symmetric
+        copyUpperToLower(4, m);
+        m *= k/6.;
+    }
+};
+
+} // namespace
+
+#endif
+
diff --git a/Tests/NumLib/TestCoordinatesMapping.cpp b/Tests/NumLib/TestCoordinatesMapping.cpp
index 4d5739b4d7ed7b505d904022b121c1197b97396e..f4a5bf67f8985269b7a3478f84992a79418810e8 100644
--- a/Tests/NumLib/TestCoordinatesMapping.cpp
+++ b/Tests/NumLib/TestCoordinatesMapping.cpp
@@ -1,7 +1,4 @@
 /**
- * \author Norihiro Watanabe
- * \date   2013-09-06
- *
  * \copyright
  * Copyright (c) 2013, OpenGeoSys Community (http://www.opengeosys.org)
  *            Distributed under a Modified BSD License.
@@ -19,53 +16,70 @@
 #include <Eigen/Eigen>
 #endif
 
-#include "MeshLib/Elements/Quad.h"
-#include "NumLib/Fem/ShapeFunction/ShapeQuad4.h"
 #include "NumLib/Fem/CoordinatesMapping/ShapeMatrices.h"
 #include "NumLib/Fem/CoordinatesMapping/NaturalCoordinatesMapping.h"
 
+#include "CoordinatesMappingTestData/TestLine2.h"
+#include "CoordinatesMappingTestData/TestQuad4.h"
+#include "CoordinatesMappingTestData/TestHex8.h"
+
 #include "../TestTools.h"
 
 using namespace NumLib;
+using namespace CoordinatesMappingTestData;
 
-#ifdef OGS_USE_EIGEN
 
 namespace
 {
-class NumLibFemNaturalCoordinatesMappingQuad4Test : public ::testing::Test
+
+// Element types to be tested
+typedef ::testing::Types<
+        TestLine2,
+        TestQuad4,
+        TestHex8
+        > TestTypes;
+}
+
+template <class T_TEST>
+class NumLibFemNaturalCoordinatesMappingTest : public ::testing::Test, public T_TEST
 {
- public:
+public:
+    typedef typename T_TEST::ElementType ElementType;
+    typedef typename T_TEST::ShapeFunctionType ShapeFunctionType;
+    static const unsigned dim = T_TEST::dim;
+    static const unsigned e_nnodes = T_TEST::e_nnodes;
     // Matrix types
-    static const unsigned dim = 2;
-    static const unsigned e_nnodes = 4;
+#ifdef OGS_USE_EIGEN
     typedef Eigen::Matrix<double, e_nnodes, 1> NodalVector;
     typedef Eigen::Matrix<double, dim, e_nnodes, Eigen::RowMajor> DimNodalMatrix;
     typedef Eigen::Matrix<double, dim, dim, Eigen::RowMajor> DimMatrix;
+#endif
     // Shape data type
     typedef ShapeMatrices<NodalVector,DimNodalMatrix,DimMatrix> ShapeMatricesType;
     // Natural coordinates mapping type
-    typedef NaturalCoordinatesMapping<MeshLib::Quad, ShapeQuad4, ShapeMatricesType> NaturalCoordsMappingType;
+    typedef NaturalCoordinatesMapping<ElementType, ShapeFunctionType, ShapeMatricesType> NaturalCoordsMappingType;
 
- public:
-    NumLibFemNaturalCoordinatesMappingQuad4Test()
+public:
+    NumLibFemNaturalCoordinatesMappingTest()
+    : eps(std::numeric_limits<double>::epsilon())
     {
-        // create four quad elements used for testing
-        naturalQuad   = createNaturalShapeQuad();
-        irregularQuad = createIrregularShapeQuad();
-        clockwiseQuad = createClockWiseQuad();
-        zeroAreaQuad  = createZeroAreaQuad();
+        // create four elements used for testing
+        naturalEle   = this->createNaturalShape();
+        irregularEle = this->createIrregularShape();
+        clockwiseEle  = this->createClockWise();
+        zeroVolumeEle  = this->createZeroVolume();
 
         // for destructor
-        vec_eles.push_back(naturalQuad);
-        vec_eles.push_back(irregularQuad);
-        vec_eles.push_back(clockwiseQuad);
-        vec_eles.push_back(zeroAreaQuad);
+        vec_eles.push_back(naturalEle);
+        vec_eles.push_back(irregularEle);
+        vec_eles.push_back(clockwiseEle);
+        vec_eles.push_back(zeroVolumeEle);
         for (auto e : vec_eles)
             for (unsigned i=0; i<e->getNNodes(true); i++)
                 vec_nodes.push_back(e->getNode(i));
     }
 
-    ~NumLibFemNaturalCoordinatesMappingQuad4Test()
+    ~NumLibFemNaturalCoordinatesMappingTest()
     {
         for (auto itr = vec_nodes.begin(); itr!=vec_nodes.end(); ++itr )
             delete *itr;
@@ -73,77 +87,25 @@ class NumLibFemNaturalCoordinatesMappingQuad4Test : public ::testing::Test
             delete *itr;
     }
 
-    // quad having shape identical to that in natural coordinates
-    MeshLib::Quad* createNaturalShapeQuad()
-    {
-        MeshLib::Node** nodes = new MeshLib::Node*[e_nnodes];
-        nodes[0] = new MeshLib::Node( 1.0,  1.0,  0.0);
-        nodes[1] = new MeshLib::Node(-1.0,  1.0,  0.0);
-        nodes[2] = new MeshLib::Node(-1.0, -1.0,  0.0);
-        nodes[3] = new MeshLib::Node( 1.0, -1.0,  0.0);
-        return new MeshLib::Quad(nodes);
-    }
-
-    // quad having irregular or skew shape
-    MeshLib::Quad* createIrregularShapeQuad()
-    {
-        MeshLib::Node** nodes = new MeshLib::Node*[e_nnodes];
-        nodes[0] = new MeshLib::Node(-0.5, -0.5,  0.0);
-        nodes[1] = new MeshLib::Node( 0.6, -0.6,  0.0);
-        nodes[2] = new MeshLib::Node( 0.5,  0.4,  0.0);
-        nodes[3] = new MeshLib::Node(-0.3,  0.1,  0.0);
-        return new MeshLib::Quad(nodes);
-    }
-
-    // invalid case: clock wise node ordering
-    MeshLib::Quad* createClockWiseQuad()
-    {
-        MeshLib::Node** nodes = new MeshLib::Node*[e_nnodes];
-        nodes[0] = new MeshLib::Node( 1.0,  1.0,  0.0);
-        nodes[3] = new MeshLib::Node(-1.0,  1.0,  0.0);
-        nodes[2] = new MeshLib::Node(-1.0, -1.0,  0.0);
-        nodes[1] = new MeshLib::Node( 1.0, -1.0,  0.0);
-        return new MeshLib::Quad(nodes);
-    }
-
-    // invalid case: zero area
-    MeshLib::Quad* createZeroAreaQuad()
-    {
-        MeshLib::Node** nodes = new MeshLib::Node*[e_nnodes];
-        nodes[0] = new MeshLib::Node( 1.0,  1.0,  0.0);
-        nodes[1] = new MeshLib::Node(-1.0,  1.0,  0.0);
-        nodes[2] = new MeshLib::Node(-1.0,  1.0,  0.0);
-        nodes[3] = new MeshLib::Node( 1.0,  1.0,  0.0);
-        return new MeshLib::Quad(nodes);
-    }
-
-    static const double r[dim];
-    static const double exp_N[e_nnodes];
-    static const double exp_dNdr[e_nnodes*dim];
-    static const double eps;
-
+    const double eps;
     std::vector<const MeshLib::Node*> vec_nodes;
-    std::vector<const MeshLib::Quad*> vec_eles;
-    MeshLib::Quad* naturalQuad;
-    MeshLib::Quad* irregularQuad;
-    MeshLib::Quad* clockwiseQuad;
-    MeshLib::Quad* zeroAreaQuad;
-
-}; // NumLibFemNaturalCoordinatesMappingQuad4Test
+    std::vector<const ElementType*> vec_eles;
+    ElementType* naturalEle;
+    ElementType* irregularEle;
+    ElementType* clockwiseEle;
+    ElementType* zeroVolumeEle;
+};
 
-const double NumLibFemNaturalCoordinatesMappingQuad4Test::r[dim] = {0.5, 0.5};
-const double NumLibFemNaturalCoordinatesMappingQuad4Test::exp_N[e_nnodes] = {0.5625, 0.1875, 0.0625, 0.1875};
-const double NumLibFemNaturalCoordinatesMappingQuad4Test::exp_dNdr[e_nnodes*dim] = {0.375, -0.375, -0.125, 0.125, 0.375, 0.125, -0.125, -0.375};
-const double NumLibFemNaturalCoordinatesMappingQuad4Test::eps = std::numeric_limits<double>::epsilon();
+TYPED_TEST_CASE(NumLibFemNaturalCoordinatesMappingTest, TestTypes);
 
-} // namespace
-
-TEST_F(NumLibFemNaturalCoordinatesMappingQuad4Test, CheckFieldSpecification_N)
+TYPED_TEST(NumLibFemNaturalCoordinatesMappingTest, CheckFieldSpecification_N)
 {
-    ShapeMatricesType shape(dim, e_nnodes);
+    typedef typename TestFixture::ShapeMatricesType ShapeMatricesType;
+    typedef typename TestFixture::NaturalCoordsMappingType NaturalCoordsMappingType;
+    ShapeMatricesType shape(this->dim, this->e_nnodes);
 
     //only N
-    NaturalCoordsMappingType::computeShapeMatrices<ShapeMatrixType::N>(*naturalQuad, r, shape);
+    NaturalCoordsMappingType::template computeShapeMatrices<ShapeMatrixType::N>(*this->naturalEle, this->r, shape);
     ASSERT_FALSE(shape.N.isZero());
     ASSERT_TRUE(shape.dNdr.isZero());
     ASSERT_TRUE(shape.J.isZero());
@@ -152,12 +114,14 @@ TEST_F(NumLibFemNaturalCoordinatesMappingQuad4Test, CheckFieldSpecification_N)
     ASSERT_TRUE(shape.dNdx.isZero());
 }
 
-TEST_F(NumLibFemNaturalCoordinatesMappingQuad4Test, CheckFieldSpecification_DNDR)
+TYPED_TEST(NumLibFemNaturalCoordinatesMappingTest, CheckFieldSpecification_DNDR)
 {
-    ShapeMatricesType shape(dim, e_nnodes);
+    typedef typename TestFixture::ShapeMatricesType ShapeMatricesType;
+    typedef typename TestFixture::NaturalCoordsMappingType NaturalCoordsMappingType;
+    ShapeMatricesType shape(this->dim, this->e_nnodes);
 
     // dNdr
-    NaturalCoordsMappingType::computeShapeMatrices<ShapeMatrixType::DNDR>(*naturalQuad, r, shape);
+    NaturalCoordsMappingType::template computeShapeMatrices<ShapeMatrixType::DNDR>(*this->naturalEle, this->r, shape);
     ASSERT_TRUE(shape.N.isZero());
     ASSERT_FALSE(shape.dNdr.isZero());
     ASSERT_TRUE(shape.J.isZero());
@@ -166,13 +130,15 @@ TEST_F(NumLibFemNaturalCoordinatesMappingQuad4Test, CheckFieldSpecification_DNDR
     ASSERT_TRUE(shape.dNdx.isZero());
 }
 
-TEST_F(NumLibFemNaturalCoordinatesMappingQuad4Test, CheckFieldSpecification_N_J)
+TYPED_TEST(NumLibFemNaturalCoordinatesMappingTest, CheckFieldSpecification_N_J)
 {
-    ShapeMatricesType shape(dim, e_nnodes);
+    typedef typename TestFixture::ShapeMatricesType ShapeMatricesType;
+    typedef typename TestFixture::NaturalCoordsMappingType NaturalCoordsMappingType;
+    ShapeMatricesType shape(this->dim, this->e_nnodes);
 
     // N_J
     shape.setZero();
-    NaturalCoordsMappingType::computeShapeMatrices<ShapeMatrixType::N_J>(*naturalQuad, r, shape);
+    NaturalCoordsMappingType::template computeShapeMatrices<ShapeMatrixType::N_J>(*this->naturalEle, this->r, shape);
     ASSERT_FALSE(shape.N.isZero());
     ASSERT_FALSE(shape.dNdr.isZero());
     ASSERT_FALSE(shape.J.isZero());
@@ -181,12 +147,14 @@ TEST_F(NumLibFemNaturalCoordinatesMappingQuad4Test, CheckFieldSpecification_N_J)
     ASSERT_TRUE(shape.dNdx.isZero());
 }
 
-TEST_F(NumLibFemNaturalCoordinatesMappingQuad4Test, CheckFieldSpecification_DNDR_J)
+TYPED_TEST(NumLibFemNaturalCoordinatesMappingTest, CheckFieldSpecification_DNDR_J)
 {
-    ShapeMatricesType shape(dim, e_nnodes);
+    typedef typename TestFixture::ShapeMatricesType ShapeMatricesType;
+    typedef typename TestFixture::NaturalCoordsMappingType NaturalCoordsMappingType;
+    ShapeMatricesType shape(this->dim, this->e_nnodes);
 
     // dNdr, J
-    NaturalCoordsMappingType::computeShapeMatrices<ShapeMatrixType::DNDR_J>(*naturalQuad, r, shape);
+    NaturalCoordsMappingType::template computeShapeMatrices<ShapeMatrixType::DNDR_J>(*this->naturalEle, this->r, shape);
     ASSERT_TRUE(shape.N.isZero());
     ASSERT_FALSE(shape.dNdr.isZero());
     ASSERT_FALSE(shape.J.isZero());
@@ -195,13 +163,15 @@ TEST_F(NumLibFemNaturalCoordinatesMappingQuad4Test, CheckFieldSpecification_DNDR
     ASSERT_TRUE(shape.dNdx.isZero());
 }
 
-TEST_F(NumLibFemNaturalCoordinatesMappingQuad4Test, CheckFieldSpecification_DNDX)
+TYPED_TEST(NumLibFemNaturalCoordinatesMappingTest, CheckFieldSpecification_DNDX)
 {
-    ShapeMatricesType shape(dim, e_nnodes);
+    typedef typename TestFixture::ShapeMatricesType ShapeMatricesType;
+    typedef typename TestFixture::NaturalCoordsMappingType NaturalCoordsMappingType;
+    ShapeMatricesType shape(this->dim, this->e_nnodes);
 
     // DNDX
     shape.setZero();
-    NaturalCoordsMappingType::computeShapeMatrices<ShapeMatrixType::DNDX>(*naturalQuad, r, shape);
+    NaturalCoordsMappingType::template computeShapeMatrices<ShapeMatrixType::DNDX>(*this->naturalEle, this->r, shape);
     ASSERT_TRUE(shape.N.isZero());
     ASSERT_FALSE(shape.dNdr.isZero());
     ASSERT_FALSE(shape.J.isZero());
@@ -210,13 +180,15 @@ TEST_F(NumLibFemNaturalCoordinatesMappingQuad4Test, CheckFieldSpecification_DNDX
     ASSERT_FALSE(shape.dNdx.isZero());
 }
 
-TEST_F(NumLibFemNaturalCoordinatesMappingQuad4Test, CheckFieldSpecification_ALL)
+TYPED_TEST(NumLibFemNaturalCoordinatesMappingTest, CheckFieldSpecification_ALL)
 {
-    ShapeMatricesType shape(dim, e_nnodes);
+    typedef typename TestFixture::ShapeMatricesType ShapeMatricesType;
+    typedef typename TestFixture::NaturalCoordsMappingType NaturalCoordsMappingType;
+    ShapeMatricesType shape(this->dim, this->e_nnodes);
 
     // ALL
     shape.setZero();
-    NaturalCoordsMappingType::computeShapeMatrices(*naturalQuad, r, shape);
+    NaturalCoordsMappingType::computeShapeMatrices(*this->naturalEle, this->r, shape);
     ASSERT_FALSE(shape.N.isZero());
     ASSERT_FALSE(shape.dNdr.isZero());
     ASSERT_FALSE(shape.J.isZero());
@@ -225,82 +197,83 @@ TEST_F(NumLibFemNaturalCoordinatesMappingQuad4Test, CheckFieldSpecification_ALL)
     ASSERT_FALSE(shape.dNdx.isZero());
 }
 
-
-TEST_F(NumLibFemNaturalCoordinatesMappingQuad4Test, CheckNaturalShape)
+TYPED_TEST(NumLibFemNaturalCoordinatesMappingTest, CheckNaturalShape)
 {
-    // identical to natural coordinates
-    ShapeMatricesType shape(dim, e_nnodes);
-
-    NaturalCoordsMappingType::computeShapeMatrices(*naturalQuad, r, shape);
-    double exp_J[]= {1.0, 0.0, 0.0, 1.0};
+    typedef typename TestFixture::ShapeMatricesType ShapeMatricesType;
+    typedef typename TestFixture::NaturalCoordsMappingType NaturalCoordsMappingType;
 
-    ASSERT_ARRAY_NEAR(exp_N, shape.N.data(), shape.N.size(), eps);
-    ASSERT_ARRAY_NEAR(exp_dNdr, shape.dNdr.data(), shape.dNdr.size(), eps);
-    ASSERT_ARRAY_NEAR(exp_J, shape.J.data(), shape.J.size(), eps);
-    ASSERT_ARRAY_NEAR(exp_J, shape.invJ.data(), shape.invJ.size(), eps);
-    ASSERT_NEAR(1.0, shape.detJ, eps);
-    ASSERT_ARRAY_NEAR(exp_dNdr, shape.dNdx.data(), shape.dNdx.size(), eps);
+    // identical to natural coordinates
+    ShapeMatricesType shape(this->dim, this->e_nnodes);
+    NaturalCoordsMappingType::computeShapeMatrices(*this->naturalEle, this->r, shape);
+    double exp_J[TestFixture::dim*TestFixture::dim]= {0.0};
+    for (unsigned i=0; i<this->dim; i++)
+        exp_J[i+this->dim*i] = 1.0;
+
+    ASSERT_ARRAY_NEAR(this->nat_exp_N, shape.N.data(), shape.N.size(), this->eps);
+    ASSERT_ARRAY_NEAR(this->nat_exp_dNdr, shape.dNdr.data(), shape.dNdr.size(), this->eps);
+    ASSERT_ARRAY_NEAR(exp_J, shape.J.data(), shape.J.size(), this->eps);
+    ASSERT_ARRAY_NEAR(exp_J, shape.invJ.data(), shape.invJ.size(), this->eps);
+    ASSERT_NEAR(1.0, shape.detJ, this->eps);
+    ASSERT_ARRAY_NEAR(this->nat_exp_dNdr, shape.dNdx.data(), shape.dNdx.size(), this->eps);
 }
 
-TEST_F(NumLibFemNaturalCoordinatesMappingQuad4Test, CheckIrregularShape)
+TYPED_TEST(NumLibFemNaturalCoordinatesMappingTest, CheckIrregularShape)
 {
+    typedef typename TestFixture::ShapeMatricesType ShapeMatricesType;
+    typedef typename TestFixture::NaturalCoordsMappingType NaturalCoordsMappingType;
+
     // irregular shape
-    ShapeMatricesType shape(dim, e_nnodes);
-
-    NaturalCoordsMappingType::computeShapeMatrices(*irregularQuad, r, shape);
-//        std::cout << shape;
-    double exp_J[]= {-0.5125, 0.0, -0.0625, -0.35};
-    double exp_invJ[]= {-1.9512195121951219, 0.0, 0.3484320557491290, -2.8571428571428572};
-    double exp_dNdx[]= {-0.73170731707317072, 0.73170731707317072, 0.243902439024390, -0.24390243902439029, -0.940766550522648, -0.48780487804878048, 0.313588850174216, 1.1149825783972125};
-
-    ASSERT_ARRAY_NEAR(exp_N, shape.N.data(), shape.N.size(), eps);
-    ASSERT_ARRAY_NEAR(exp_dNdr, shape.dNdr.data(), shape.dNdr.size(), eps);
-    ASSERT_ARRAY_NEAR(exp_J, shape.J.data(), shape.J.size(), eps);
-    ASSERT_ARRAY_NEAR(exp_invJ, shape.invJ.data(), shape.invJ.size(), eps);
-    ASSERT_NEAR(0.179375, shape.detJ, eps);
-    ASSERT_ARRAY_NEAR(exp_dNdx, shape.dNdx.data(), shape.dNdx.size(), eps);
+    ShapeMatricesType shape(this->dim, this->e_nnodes);
+    NaturalCoordsMappingType::computeShapeMatrices(*this->irregularEle, this->r, shape);
+    //std::cout <<  std::setprecision(16) << shape;
+
+    ASSERT_ARRAY_NEAR(this->nat_exp_N, shape.N.data(), shape.N.size(), this->eps);
+    ASSERT_ARRAY_NEAR(this->nat_exp_dNdr, shape.dNdr.data(), shape.dNdr.size(), this->eps);
+    ASSERT_ARRAY_NEAR(this->ir_exp_J, shape.J.data(), shape.J.size(), this->eps);
+    ASSERT_NEAR(this->ir_exp_detJ, shape.detJ, this->eps);
+    ASSERT_ARRAY_NEAR(this->ir_exp_invJ, shape.invJ.data(), shape.invJ.size(), this->eps);
+    ASSERT_ARRAY_NEAR(this->ir_exp_dNdx, shape.dNdx.data(), shape.dNdx.size(), this->eps);
 }
 
-TEST_F(NumLibFemNaturalCoordinatesMappingQuad4Test, CheckClockwise)
+TYPED_TEST(NumLibFemNaturalCoordinatesMappingTest, CheckClockwise)
 {
-    // clockwise node ordering, which is invalid)
-    ShapeMatricesType shape(dim, e_nnodes);
+    typedef typename TestFixture::ShapeMatricesType ShapeMatricesType;
+    typedef typename TestFixture::NaturalCoordsMappingType NaturalCoordsMappingType;
 
-    NaturalCoordsMappingType::computeShapeMatrices(*clockwiseQuad, r, shape);
-    //std::cout << shape;
-    double exp_J[]= {0.0, 1.0, 1.0, 0.0};
+    // clockwise node ordering, which is invalid)
+    ShapeMatricesType shape(this->dim, this->e_nnodes);
+    NaturalCoordsMappingType::computeShapeMatrices(*this->clockwiseEle, this->r, shape);
+    //std::cout <<  std::setprecision(16) << shape;
     // Inverse of the Jacobian matrix doesn't exist
-    double exp_invJ[dim*dim]= {0.0};
-    double exp_dNdx[dim*e_nnodes]= {0.0};
-
-    ASSERT_ARRAY_NEAR(exp_N, shape.N.data(), shape.N.size(), eps);
-    ASSERT_ARRAY_NEAR(exp_dNdr, shape.dNdr.data(), shape.dNdr.size(), eps);
-    ASSERT_ARRAY_NEAR(exp_J, shape.J.data(), shape.J.size(), eps);
-    ASSERT_ARRAY_NEAR(exp_invJ, shape.invJ.data(), shape.invJ.size(), eps);
-    ASSERT_NEAR(-1.0, shape.detJ, eps);
-    ASSERT_ARRAY_NEAR(exp_dNdx, shape.dNdx.data(), shape.dNdx.size(), eps);
+    double exp_invJ[TestFixture::dim*TestFixture::dim]= {0.0};
+    double exp_dNdx[TestFixture::dim*TestFixture::e_nnodes]= {0.0};
+
+    ASSERT_ARRAY_NEAR(this->nat_exp_N, shape.N.data(), shape.N.size(), this->eps);
+    ASSERT_ARRAY_NEAR(this->nat_exp_dNdr, shape.dNdr.data(), shape.dNdr.size(), this->eps);
+    ASSERT_ARRAY_NEAR(this->cl_exp_J, shape.J.data(), shape.J.size(), this->eps);
+    ASSERT_NEAR(this->cl_exp_detJ, shape.detJ, this->eps);
+    ASSERT_ARRAY_NEAR(exp_invJ, shape.invJ.data(), shape.invJ.size(), this->eps);
+    ASSERT_ARRAY_NEAR(exp_dNdx, shape.dNdx.data(), shape.dNdx.size(), this->eps);
 }
 
-TEST_F(NumLibFemNaturalCoordinatesMappingQuad4Test, CheckZeroArea)
+TYPED_TEST(NumLibFemNaturalCoordinatesMappingTest, CheckZeroVolume)
 {
-    // zero area
-    ShapeMatricesType shape(dim, e_nnodes);
+    typedef typename TestFixture::ShapeMatricesType ShapeMatricesType;
+    typedef typename TestFixture::NaturalCoordsMappingType NaturalCoordsMappingType;
 
-    NaturalCoordsMappingType::computeShapeMatrices(*zeroAreaQuad, r, shape);
-    //std::cout << shape;
-    double exp_J[]= {1.0, 0.0, 0.0, 0.0};
+    ShapeMatricesType shape(this->dim, this->e_nnodes);
+    NaturalCoordsMappingType::computeShapeMatrices(*this->zeroVolumeEle, this->r, shape);
+    //std::cout <<  std::setprecision(16) << shape;
     // Inverse of the Jacobian matrix doesn't exist
-    double exp_invJ[dim*dim]= {0.0};
-    double exp_dNdx[dim*e_nnodes]= {0.0};
-
-    ASSERT_ARRAY_NEAR(exp_N, shape.N.data(), shape.N.size(), eps);
-    ASSERT_ARRAY_NEAR(exp_dNdr, shape.dNdr.data(), shape.dNdr.size(), eps);
-    ASSERT_ARRAY_NEAR(exp_J, shape.J.data(), shape.J.size(), eps);
-    ASSERT_ARRAY_NEAR(exp_invJ, shape.invJ.data(), shape.invJ.size(), eps);
-    ASSERT_NEAR(0.0, shape.detJ, eps);
-    ASSERT_ARRAY_NEAR(exp_dNdx, shape.dNdx.data(), shape.dNdx.size(), eps);
+    double exp_invJ[TestFixture::dim*TestFixture::dim]= {0.0};
+    double exp_dNdx[TestFixture::dim*TestFixture::e_nnodes]= {0.0};
+
+    ASSERT_ARRAY_NEAR(this->nat_exp_N, shape.N.data(), shape.N.size(), this->eps);
+    ASSERT_ARRAY_NEAR(this->nat_exp_dNdr, shape.dNdr.data(), shape.dNdr.size(), this->eps);
+    ASSERT_ARRAY_NEAR(this->ze_exp_J, shape.J.data(), shape.J.size(), this->eps);
+    ASSERT_NEAR(0.0, shape.detJ, this->eps);
+    ASSERT_ARRAY_NEAR(exp_invJ, shape.invJ.data(), shape.invJ.size(), this->eps);
+    ASSERT_ARRAY_NEAR(exp_dNdx, shape.dNdx.data(), shape.dNdx.size(), this->eps);
 }
 
-#endif //OGS_USE_EIGEN
-
 
diff --git a/Tests/NumLib/TestFeQuad4.cpp b/Tests/NumLib/TestFe.cpp
similarity index 54%
rename from Tests/NumLib/TestFeQuad4.cpp
rename to Tests/NumLib/TestFe.cpp
index 8f47b0b7f78cd706252ffb14c71dc1804d1ecb87..13b9fd48cf976f0ac91d7eb5d61b951184f86960 100644
--- a/Tests/NumLib/TestFeQuad4.cpp
+++ b/Tests/NumLib/TestFe.cpp
@@ -1,8 +1,4 @@
 /**
- * \file TestFEM.cpp
- * \author Norihiro Watanabe
- * \date   2012-08-03
- *
  * \copyright
  * Copyright (c) 2013, OpenGeoSys Community (http://www.opengeosys.org)
  *            Distributed under a Modified BSD License.
@@ -19,60 +15,110 @@
 #include <Eigen/Eigen>
 #endif
 
-#include "MeshLib/Elements/Quad.h"
 #include "NumLib/Fem/CoordinatesMapping/ShapeMatrices.h"
 #include "NumLib/Fem/FiniteElement/C0IsoparametricElements.h"
 
+#include "FeTestData/TestFeLINE2.h"
+#include "FeTestData/TestFeQUAD4.h"
+#include "FeTestData/TestFeHEX8.h"
+
 #include "Tests/TestTools.h"
 
 using namespace NumLib;
+using namespace FeTestData;
 
 namespace
 {
+#ifdef OGS_USE_EIGEN
+template <typename T_FE>
+struct EigenFixedMatrixTypes
+{
+    typedef Eigen::Matrix<double, T_FE::e_nnodes, T_FE::e_nnodes, Eigen::RowMajor> NodalMatrixType;
+    typedef Eigen::Matrix<double, T_FE::e_nnodes, 1> NodalVectorType;
+    typedef Eigen::Matrix<double, T_FE::dim, T_FE::e_nnodes, Eigen::RowMajor> DimNodalMatrixType;
+    typedef Eigen::Matrix<double, T_FE::dim, T_FE::dim, Eigen::RowMajor> DimMatrixType;
+};
+
+struct EigenDynamicMatrixTypes
+{
+    typedef Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor> NodalMatrixType;
+    typedef NodalMatrixType DimNodalMatrixType;
+    typedef NodalMatrixType DimMatrixType;
+    typedef Eigen::VectorXd NodalVectorType;
+};
+
+typedef EigenDynamicMatrixTypes DefaultMatrixType;
+#endif // OGS_USE_EIGEN
+
+// test cases
+template <class T_FE_TYPE, class T_MAT_TYPES=DefaultMatrixType>
+struct TestCase
+{
+    typedef T_FE_TYPE T_FE;
+    typedef T_MAT_TYPES T_MATRIX_TYPES;
+};
 
-static const unsigned dim = 2;
-static const unsigned e_nnodes = 4;
+typedef ::testing::Types<
+        TestCase<TestFeLINE2>,
+        TestCase<TestFeQUAD4>,
+        TestCase<TestFeHEX8>
+#ifdef OGS_USE_EIGEN
+        ,TestCase<TestFeLINE2, EigenFixedMatrixTypes<TestFeLINE2> >
+        ,TestCase<TestFeQUAD4, EigenFixedMatrixTypes<TestFeQUAD4> >
+        ,TestCase<TestFeHEX8, EigenFixedMatrixTypes<TestFeHEX8> >
+#endif
+        > TestTypes;
+}
 
-template <class T_MATRIX_TYPES>
-class NumLibFemIsoQuad4Test : public ::testing::Test
+template <class T>
+class NumLibFemIsoTest : public ::testing::Test, public T::T_FE
 {
  public:
+    typedef typename T::T_MATRIX_TYPES T_MATRIX_TYPES;
+    typedef typename T::T_FE T_FE;
     // Matrix types
+    typedef T_MATRIX_TYPES MatrixType;
     typedef typename T_MATRIX_TYPES::NodalMatrixType NodalMatrix;
     typedef typename T_MATRIX_TYPES::NodalVectorType NodalVector;
     typedef typename T_MATRIX_TYPES::DimNodalMatrixType DimNodalMatrix;
     typedef typename T_MATRIX_TYPES::DimMatrixType DimMatrix;
     // Finite element type
-    typedef typename NumLib::FeQUAD4<NodalVector, DimNodalMatrix, DimMatrix>::type FeQUAD4Type;
+    typedef typename T_FE::template FeType<T_MATRIX_TYPES>::type::type FeType;
     // Shape matrix data type
-    typedef typename FeQUAD4Type::ShapeMatricesType ShapeMatricesType;
+    typedef typename FeType::ShapeMatricesType ShapeMatricesType;
+    typedef typename T_FE::MeshElementType MeshElementType;
+
+    static const unsigned dim = T_FE::dim;
+    static const unsigned e_nnodes = T_FE::e_nnodes;
+    static const unsigned n_sample_pt_order2 = T_FE::n_sample_pt_order2;
+    static const unsigned n_sample_pt_order3 = T_FE::n_sample_pt_order3;
 
  public:
-    NumLibFemIsoQuad4Test() :
+    NumLibFemIsoTest() :
         D(dim, dim),
         expectedM(e_nnodes,e_nnodes),
         expectedK(e_nnodes,e_nnodes),
         integration_method(2)
     {
-        // create a quad element used for testing
-        unitSquareQuad   = createSquareQuad(1.0);
+        // create a mesh element used for testing
+        mesh_element = this->createMeshElement();
 
         // set a conductivity tensor
         setIdentityMatrix(dim, D);
         D *= conductivity;
 
         // set expected matrices
-        setExpectedMassMatrix(expectedM);
-        setExpectedLaplaceMatrix(conductivity, expectedK);
+        this->setExpectedMassMatrix(expectedM);
+        this->setExpectedLaplaceMatrix(conductivity, expectedK);
 
         // for destructor
-        vec_eles.push_back(unitSquareQuad);
+        vec_eles.push_back(mesh_element);
         for (auto e : vec_eles)
             for (unsigned i=0; i<e->getNNodes(true); i++)
                 vec_nodes.push_back(e->getNode(i));
     }
 
-    ~NumLibFemIsoQuad4Test()
+    virtual ~NumLibFemIsoTest()
     {
         for (auto itr = vec_nodes.begin(); itr!=vec_nodes.end(); ++itr )
             delete *itr;
@@ -80,128 +126,54 @@ class NumLibFemIsoQuad4Test : public ::testing::Test
             delete *itr;
     }
 
-    // Quad: square shape with a length of 1m
-    MeshLib::Quad* createSquareQuad(double h)
-    {
-        MeshLib::Node** nodes = new MeshLib::Node*[e_nnodes];
-        nodes[0] = new MeshLib::Node(0.0, 0.0, 0.0);
-        nodes[1] = new MeshLib::Node(  h, 0.0, 0.0);
-        nodes[2] = new MeshLib::Node(  h,   h, 0.0);
-        nodes[3] = new MeshLib::Node(0.0,   h, 0.0);
-        return new MeshLib::Quad(nodes);
-    }
-
-    // set an identity matrix
-    template <class T_MATRIX, typename ID_TYPE=signed>
-    void setIdentityMatrix(unsigned dim, T_MATRIX &m) const
-    {
-        for (unsigned i=0; i<dim; i++)
-            for (unsigned j=0; j<dim; j++)
-                m(i,j) = 0.0;
-        for (unsigned i=0; i<dim; i++)
-            m(i,i) = 1.0;
-    }
-
-    // copy upper triangles to lower triangles in a matrix
-    template <class T_MATRIX, typename ID_TYPE=signed>
-    void copyUpperToLower(const ID_TYPE dim, T_MATRIX &m) const
-    {
-        for (ID_TYPE i=0; i<dim; i++)
-            for (ID_TYPE j=0; j<i; j++)
-                m(i,j) = m(j,i);
-    }
-
-    // set an expected mass matrix for 1m x 1m
-    template <class T_MATRIX, typename ID_TYPE=signed>
-    void setExpectedMassMatrix(T_MATRIX &m)
-    {
-        // set upper triangle entries
-        m(0,0) = 1.0; m(0,1) = 1./2; m(0,2) = 1./4; m(0,3) = 1./2;
-        m(1,1) = 1.0; m(1,2) = 1./2; m(1,3) = 1./4;
-        m(2,2) = 1.0; m(2,3) = 1./2;
-        m(3,3) = 1.0;
-        // make symmetric
-        copyUpperToLower(4, m);
-        m *= 1./9.;
-    }
-
-    // set an expected laplace matrix for 1m x 1m
-    template <class T_MATRIX, typename ID_TYPE=signed>
-    void setExpectedLaplaceMatrix(double k, T_MATRIX &m)
-    {
-        // set upper triangle entries
-        m(0,0) = 4.0; m(0,1) = -1.0; m(0,2) = -2.0; m(0,3) = -1.0;
-        m(1,1) = 4.0; m(1,2) = -1.0; m(1,3) = -2.0;
-        m(2,2) = 4.0; m(2,3) = -1.0;
-        m(3,3) = 4.0;
-        // make symmetric
-        copyUpperToLower(4, m);
-        m *= k/6.;
-    }
 
     static const double conductivity;
     static const double eps;
     DimMatrix D;
     NodalMatrix expectedM;
     NodalMatrix expectedK;
-    typename FeQUAD4Type::IntegrationMethod integration_method;
+    typename FeType::IntegrationMethod integration_method;
 
     std::vector<const MeshLib::Node*> vec_nodes;
-    std::vector<const MeshLib::Quad*> vec_eles;
-    MeshLib::Quad* unitSquareQuad;
+    std::vector<const MeshElementType*> vec_eles;
+    MeshElementType* mesh_element;
 
-}; // NumLibFemIsoQuad4Test
+}; // NumLibFemIsoTest
 
-template <class T_MATRIX_TYPES>
-const double NumLibFemIsoQuad4Test<T_MATRIX_TYPES>::conductivity = 1e-11;
+template <class T>
+const double NumLibFemIsoTest<T>::conductivity = 1e-11;
 
-template <class T_MATRIX_TYPES>
-const double NumLibFemIsoQuad4Test<T_MATRIX_TYPES>::eps = std::numeric_limits<double>::epsilon();
+template <class T>
+const double NumLibFemIsoTest<T>::eps = std::numeric_limits<double>::epsilon();
 
-} // namespace
+template <class T>
+const unsigned NumLibFemIsoTest<T>::dim;
 
-#ifdef OGS_USE_EIGEN
+template <class T>
+const unsigned NumLibFemIsoTest<T>::e_nnodes;
 
-struct EigenFixedMatrixTypes
-{
-    typedef Eigen::Matrix<double, e_nnodes, e_nnodes, Eigen::RowMajor> NodalMatrixType;
-    typedef Eigen::Matrix<double, e_nnodes, 1> NodalVectorType;
-    typedef Eigen::Matrix<double, dim, e_nnodes, Eigen::RowMajor> DimNodalMatrixType;
-    typedef Eigen::Matrix<double, dim, dim, Eigen::RowMajor> DimMatrixType;
-};
+template <class T>
+const unsigned NumLibFemIsoTest<T>::n_sample_pt_order2;
 
-struct EigenDynamicMatrixTypes
-{
-    typedef Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor> NodalMatrixType;
-    typedef NodalMatrixType DimNodalMatrixType;
-    typedef NodalMatrixType DimMatrixType;
-    typedef Eigen::VectorXd NodalVectorType;
-};
-
-#endif // OGS_USE_EIGEN
-
-typedef ::testing::Types<
-#ifdef OGS_USE_EIGEN
-        EigenFixedMatrixTypes
-        , EigenDynamicMatrixTypes
-#endif
-        > MatrixTypes;
+template <class T>
+const unsigned NumLibFemIsoTest<T>::n_sample_pt_order3;
 
-TYPED_TEST_CASE(NumLibFemIsoQuad4Test, MatrixTypes);
+TYPED_TEST_CASE(NumLibFemIsoTest, TestTypes);
 
-TYPED_TEST(NumLibFemIsoQuad4Test, CheckMassMatrix)
+TYPED_TEST(NumLibFemIsoTest, CheckMassMatrix)
 {
     // Refer to typedefs in the fixture
-    typedef typename TestFixture::FeQUAD4Type FeQUAD4Type;
+    typedef typename TestFixture::FeType FeType;
     typedef typename TestFixture::NodalMatrix NodalMatrix;
     typedef typename TestFixture::ShapeMatricesType ShapeMatricesType;
 
     // create a finite element object
-    FeQUAD4Type fe(*this->unitSquareQuad);
+    FeType fe(*this->mesh_element);
 
     // evaluate a mass matrix M = int{ N^T D N }dA_e
-    NodalMatrix M(e_nnodes, e_nnodes);
-    ShapeMatricesType shape(dim, e_nnodes);
+    NodalMatrix M(this->e_nnodes, this->e_nnodes);
+    M.setZero();
+    ShapeMatricesType shape(this->dim, this->e_nnodes);
     for (std::size_t i=0; i < this->integration_method.getNPoints(); i++) {
         shape.setZero();
         auto wp = this->integration_method.getWeightedPoint(i);
@@ -212,19 +184,20 @@ TYPED_TEST(NumLibFemIsoQuad4Test, CheckMassMatrix)
     ASSERT_ARRAY_NEAR(this->expectedM.data(), M.data(), M.size(), this->eps);
 }
 
-TYPED_TEST(NumLibFemIsoQuad4Test, CheckLaplaceMatrix)
+TYPED_TEST(NumLibFemIsoTest, CheckLaplaceMatrix)
 {
     // Refer to typedefs in the fixture
-    typedef typename TestFixture::FeQUAD4Type FeQUAD4Type;
+    typedef typename TestFixture::FeType FeType;
     typedef typename TestFixture::NodalMatrix NodalMatrix;
     typedef typename TestFixture::ShapeMatricesType ShapeMatricesType;
 
     // create a finite element object
-    FeQUAD4Type fe(*this->unitSquareQuad);
+    FeType fe(*this->mesh_element);
 
     // evaluate a Laplace matrix K = int{ dNdx^T D dNdx }dA_e
-    NodalMatrix K(e_nnodes, e_nnodes);
-    ShapeMatricesType shape(dim, e_nnodes);
+    NodalMatrix K(this->e_nnodes, this->e_nnodes);
+    K.setZero();
+    ShapeMatricesType shape(this->dim, this->e_nnodes);
     for (std::size_t i=0; i < this->integration_method.getNPoints(); i++) {
         shape.setZero();
         auto wp = this->integration_method.getWeightedPoint(i);
@@ -234,20 +207,22 @@ TYPED_TEST(NumLibFemIsoQuad4Test, CheckLaplaceMatrix)
     ASSERT_ARRAY_NEAR(this->expectedK.data(), K.data(), K.size(), this->eps);
 }
 
-TYPED_TEST(NumLibFemIsoQuad4Test, CheckMassLaplaceMatrices)
+TYPED_TEST(NumLibFemIsoTest, CheckMassLaplaceMatrices)
 {
     // Refer to typedefs in the fixture
-    typedef typename TestFixture::FeQUAD4Type FeQUAD4Type;
+    typedef typename TestFixture::FeType FeType;
     typedef typename TestFixture::NodalMatrix NodalMatrix;
     typedef typename TestFixture::ShapeMatricesType ShapeMatricesType;
 
     // create a finite element object
-    FeQUAD4Type fe(*this->unitSquareQuad);
+    FeType fe(*this->mesh_element);
 
     // evaluate both mass and laplace matrices at once
-    NodalMatrix M(e_nnodes, e_nnodes);
-    NodalMatrix K(e_nnodes, e_nnodes);
-    ShapeMatricesType shape(dim, e_nnodes);
+    NodalMatrix M(this->e_nnodes, this->e_nnodes);
+    M.setZero();
+    NodalMatrix K(this->e_nnodes, this->e_nnodes);
+    K.setZero();
+    ShapeMatricesType shape(this->dim, this->e_nnodes);
     for (std::size_t i=0; i < this->integration_method.getNPoints(); i++) {
         shape.setZero();
         auto wp = this->integration_method.getWeightedPoint(i);
@@ -259,20 +234,21 @@ TYPED_TEST(NumLibFemIsoQuad4Test, CheckMassLaplaceMatrices)
     ASSERT_ARRAY_NEAR(this->expectedK.data(), K.data(), K.size(), this->eps);
 }
 
-TYPED_TEST(NumLibFemIsoQuad4Test, CheckGaussIntegrationLevel)
+TYPED_TEST(NumLibFemIsoTest, CheckGaussIntegrationLevel)
 {
     // Refer to typedefs in the fixture
-    typedef typename TestFixture::FeQUAD4Type FeQUAD4Type;
+    typedef typename TestFixture::FeType FeType;
     typedef typename TestFixture::NodalMatrix NodalMatrix;
     typedef typename TestFixture::ShapeMatricesType ShapeMatricesType;
 
     // create a finite element object with gauss quadrature level 2
-    FeQUAD4Type fe(*this->unitSquareQuad);
+    FeType fe(*this->mesh_element);
 
     // evaluate a mass matrix
-    NodalMatrix M(e_nnodes, e_nnodes);
-    ShapeMatricesType shape(dim, e_nnodes);
-    ASSERT_EQ(4u, this->integration_method.getNPoints());
+    NodalMatrix M(this->e_nnodes, this->e_nnodes);
+    M.setZero();
+    ShapeMatricesType shape(this->dim, this->e_nnodes);
+    ASSERT_EQ(TestFixture::n_sample_pt_order2, this->integration_method.getNPoints());
     for (std::size_t i=0; i < this->integration_method.getNPoints(); i++) {
         shape.setZero();
         auto wp = this->integration_method.getWeightedPoint(i);
@@ -284,7 +260,7 @@ TYPED_TEST(NumLibFemIsoQuad4Test, CheckGaussIntegrationLevel)
     // Change gauss quadrature level to 3
     this->integration_method.setIntegrationOrder(3);
     M *= .0;
-    ASSERT_EQ(9u, this->integration_method.getNPoints());
+    ASSERT_EQ(TestFixture::n_sample_pt_order3, this->integration_method.getNPoints());
     for (std::size_t i=0; i < this->integration_method.getNPoints(); i++) {
         shape.setZero();
         auto wp = this->integration_method.getWeightedPoint(i);
diff --git a/Utils/MeshEdit/removeMeshElements.cpp b/Utils/MeshEdit/removeMeshElements.cpp
index fd33d23e5e2c570b3040670f89e62a217cf780fe..1fb36cf305df23c9e879d0f73ec65f4e3d78c84f 100644
--- a/Utils/MeshEdit/removeMeshElements.cpp
+++ b/Utils/MeshEdit/removeMeshElements.cpp
@@ -29,86 +29,7 @@
 #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;
-	}
-}
+#include "MeshEditing/ElementExtraction.h"
 
 int main (int argc, char* argv[])
 {
@@ -136,45 +57,34 @@ int main (int argc, char* argv[])
 	cmd.add(matIDArg);
 	cmd.parse(argc, argv);
 
-	MeshLib::Mesh* mesh (FileIO::readMeshFromFile(mesh_in.getValue()));
+	MeshLib::Mesh const*const mesh (FileIO::readMeshFromFile(mesh_in.getValue()));
 	INFO("Mesh read: %d nodes, %d elements.", mesh->getNNodes(), mesh->getNElements());
+	MeshLib::ElementExtraction ex(*mesh);
 
 	// 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());
+		const std::size_t n_removed_elements = ex.searchByZeroContent();
+		INFO("%d zero volume elements found.", n_removed_elements);
 	}
 	if (eleTypeArg.isSet()) {
-		std::vector<std::string> eleTypeNames = eleTypeArg.getValue();
+		const std::vector<std::string> eleTypeNames = eleTypeArg.getValue();
 		for (auto typeName : eleTypeNames) {
-			MeshElemType type = String2MeshElemType(typeName);
+			const 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());
+			const std::size_t n_removed_elements = ex.searchByElementType(type);
+			INFO("%d %s elements found.", n_removed_elements, typeName.c_str());
 		}
 	}
 	if (matIDArg.isSet()) {
-		std::vector<unsigned> vec_matID = matIDArg.getValue();
+		const 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);
+			const std::size_t n_removed_elements = ex.searchByMaterialID(matID);
+			INFO("%d elements with material ID %d found.", n_removed_elements, 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));
+	// remove the elements and create a new mesh object.
+	MeshLib::Mesh const*const new_mesh = ex.removeMeshElements(mesh->getName());
 
 	// write into a file
 	FileIO::Legacy::MeshIO meshIO;
diff --git a/Utils/SimpleMeshCreation/createMeshElemPropertiesFromASCRaster.cpp b/Utils/SimpleMeshCreation/createMeshElemPropertiesFromASCRaster.cpp
index 543d6e49dd66ce15f6608c8436d087d951e68572..60e21525704ed634514b74cc3beabfbc87e48637 100644
--- a/Utils/SimpleMeshCreation/createMeshElemPropertiesFromASCRaster.cpp
+++ b/Utils/SimpleMeshCreation/createMeshElemPropertiesFromASCRaster.cpp
@@ -142,7 +142,7 @@ int main (int argc, char* argv[])
 	}
 
 	{
-		const double mu(std::accumulate(src_properties.begin(), src_properties.end(), 0) / size);
+		const double mu(std::accumulate(src_properties.begin(), src_properties.end(), 0.0) / size);
 		INFO("Mean value of source: %f.", mu);
 
 		double src_variance(MathLib::fastpow(src_properties[0] - mu, 2));