diff --git a/Applications/CLI/Tests.cmake b/Applications/CLI/Tests.cmake
index 269adb226b6dd1d1d1309b085d09bea840ced87e..f061bbfd946a22e2211e2332d7c8a86719f67dff 100644
--- a/Applications/CLI/Tests.cmake
+++ b/Applications/CLI/Tests.cmake
@@ -99,3 +99,28 @@ foreach(mesh_size 1e5 1e6)
 		DATA square_${mesh_size}_neumann.prj square_1x1_quad_${mesh_size}.vtu square_1x1.gml
 	)
 endforeach()
+
+# LINE 1 GROUNDWATER FLOW TESTS
+foreach(mesh_size 1e1)
+	AddTest(
+		NAME GroundWaterFlowProcess_line_1_${mesh_size}
+		PATH Elliptic/line_1_GroundWaterFlow
+		EXECUTABLE ogs
+		EXECUTABLE_ARGS line_${mesh_size}.prj
+		WRAPPER time
+		TESTER vtkdiff
+		DIFF_DATA line_${mesh_size}_pcs_0_ts_1.vtu Linear_1_to_minus1 Result
+		DATA line_${mesh_size}.prj line_1_line_${mesh_size}.vtu line_1.gml
+	)
+
+	AddTest(
+        NAME GroundWaterFlowProcess_line_1_Neumann_${mesh_size}
+		PATH Elliptic/line_1_GroundWaterFlow
+		EXECUTABLE ogs
+		EXECUTABLE_ARGS line_${mesh_size}_neumann.prj
+		WRAPPER time
+		TESTER vtkdiff
+		DIFF_DATA line_${mesh_size}_neumann_pcs_0_ts_1.vtu D1_left_N1_right Result
+		DATA line_${mesh_size}_neumann.prj line_1_line_${mesh_size}.vtu line_1.gml
+	)
+endforeach()
diff --git a/AssemblerLib/LocalDataInitializer.h b/AssemblerLib/LocalDataInitializer.h
index 968c62d401ec2bb97cf47cad770a914316de2846..23a67834c04613bf705ae785f5c5399cf05f994a 100644
--- a/AssemblerLib/LocalDataInitializer.h
+++ b/AssemblerLib/LocalDataInitializer.h
@@ -23,6 +23,7 @@
 #include "NumLib/Fem/ShapeFunction/ShapeHex8.h"
 #include "NumLib/Fem/ShapeFunction/ShapeLine2.h"
 #include "NumLib/Fem/ShapeFunction/ShapeLine3.h"
+#include "NumLib/Fem/ShapeFunction/ShapePoint1.h"
 #include "NumLib/Fem/ShapeFunction/ShapePrism15.h"
 #include "NumLib/Fem/ShapeFunction/ShapePrism6.h"
 #include "NumLib/Fem/ShapeFunction/ShapePyra13.h"
@@ -73,6 +74,8 @@ public:
             [](){ return new LAData<NumLib::ShapeLine2>; };
         _builder[std::type_index(typeid(MeshLib::Line3))] =
             [](){ return new LAData<NumLib::ShapeLine3>; };
+        _builder[std::type_index(typeid(MeshLib::Point))] =
+            [](){ return new LAData<NumLib::ShapePoint1>; };
         _builder[std::type_index(typeid(MeshLib::Prism15))] =
             [](){ return new LAData<NumLib::ShapePrism15>; };
         _builder[std::type_index(typeid(MeshLib::Prism))] =
diff --git a/MeshGeoToolsLib/BoundaryElementsAtPoint.cpp b/MeshGeoToolsLib/BoundaryElementsAtPoint.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..37809ccd550b58c905988dd7dc24546340e8d243
--- /dev/null
+++ b/MeshGeoToolsLib/BoundaryElementsAtPoint.cpp
@@ -0,0 +1,41 @@
+/**
+ * @copyright
+ * Copyright (c) 2012-2015, OpenGeoSys Community (http://www.opengeosys.org)
+ *            Distributed under a Modified BSD License.
+ *              See accompanying file LICENSE.txt or
+ *              http://www.opengeosys.org/LICENSE.txt
+ */
+
+#include "BoundaryElementsAtPoint.h"
+
+#include "GeoLib/Point.h"
+#include "MeshLib/Elements/Element.h"
+#include "MeshLib/Elements/Point.h"
+#include "MeshLib/Mesh.h"
+#include "MeshLib/MeshSearch/ElementSearch.h"
+#include "MeshLib/Node.h"
+
+#include "MeshGeoToolsLib/MeshNodeSearcher.h"
+
+namespace MeshGeoToolsLib
+{
+BoundaryElementsAtPoint::BoundaryElementsAtPoint(
+    MeshLib::Mesh const &mesh, MeshNodeSearcher &mshNodeSearcher,
+    GeoLib::Point const &point)
+    : _mesh(mesh), _point(point)
+{
+	auto const node_ids = mshNodeSearcher.getMeshNodeIDs(_point);
+	assert(node_ids.size() == 1);
+	std::array<MeshLib::Node*, 1> const nodes = {{
+	    const_cast<MeshLib::Node*>(_mesh.getNode(node_ids[0]))}};
+
+	_boundary_elements.push_back(new MeshLib::Point{nodes, 0, node_ids[0]});
+}
+
+BoundaryElementsAtPoint::~BoundaryElementsAtPoint()
+{
+	for (auto p : _boundary_elements)
+		delete p;
+}
+
+}  // MeshGeoToolsLib
diff --git a/MeshGeoToolsLib/BoundaryElementsAtPoint.h b/MeshGeoToolsLib/BoundaryElementsAtPoint.h
new file mode 100644
index 0000000000000000000000000000000000000000..0a407278fd5cf644afab169fcad08fad749bfd1d
--- /dev/null
+++ b/MeshGeoToolsLib/BoundaryElementsAtPoint.h
@@ -0,0 +1,68 @@
+/**
+ * @copyright
+ * Copyright (c) 2012-2015, 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 BOUNDARYELEMENTSATPOINT_H_
+#define BOUNDARYELEMENTSATPOINT_H_
+
+#include <vector>
+
+namespace GeoLib
+{
+class Point;
+}
+
+namespace MeshLib
+{
+class Mesh;
+class Element;
+}
+
+namespace MeshGeoToolsLib
+{
+class MeshNodeSearcher;
+
+/// This class collects point elements located at a given point elements.
+class BoundaryElementsAtPoint final
+{
+public:
+	/// Constructor
+	/// \param mesh             a mesh object
+	/// \param mshNodeSearcher  a MeshNodeSearcher object which is internally
+	/// used to search mesh nodes
+	/// \param point            a point object where edges are searched
+	BoundaryElementsAtPoint(MeshLib::Mesh const& mesh,
+	                       MeshNodeSearcher& mshNodeSearcher,
+	                       GeoLib::Point const& point);
+
+	~BoundaryElementsAtPoint();
+
+	MeshLib::Mesh const& getMesh() const
+	{
+		return _mesh;
+	}
+
+	GeoLib::Point const& getPoint() const
+	{
+		return _point;
+	}
+
+	/// Return the vector of boundary elements (i.e. points).
+	std::vector<MeshLib::Element*> const& getBoundaryElements() const
+	{
+		return _boundary_elements;
+	}
+
+private:
+	MeshLib::Mesh const& _mesh;
+	GeoLib::Point const& _point;
+	std::vector<MeshLib::Element*> _boundary_elements;
+};
+
+}  // end namespace MeshGeoToolsLib
+
+#endif  // BOUNDARYELEMENTSATPOINT_H_
diff --git a/MeshGeoToolsLib/BoundaryElementsSearcher.cpp b/MeshGeoToolsLib/BoundaryElementsSearcher.cpp
index 85ef0201c60289bc51cb6727089687309429fdf3..ed168abab5dbb8c0a4c8d7634d77ab13db2f2837 100644
--- a/MeshGeoToolsLib/BoundaryElementsSearcher.cpp
+++ b/MeshGeoToolsLib/BoundaryElementsSearcher.cpp
@@ -13,9 +13,12 @@
 #include "GeoLib/Surface.h"
 
 #include "MeshLib/Mesh.h"
+#include "MeshLib/Node.h"
 #include "MeshLib/Elements/Element.h"
+#include "MeshLib/Elements/Point.h"
 
 #include "MeshGeoToolsLib/MeshNodeSearcher.h"
+#include "MeshGeoToolsLib/BoundaryElementsAtPoint.h"
 #include "MeshGeoToolsLib/BoundaryElementsAlongPolyline.h"
 #include "MeshGeoToolsLib/BoundaryElementsOnSurface.h"
 
@@ -37,6 +40,9 @@ BoundaryElementsSearcher::~BoundaryElementsSearcher()
 std::vector<MeshLib::Element*> const& BoundaryElementsSearcher::getBoundaryElements(GeoLib::GeoObject const& geoObj)
 {
 	switch (geoObj.getGeoType()) {
+	case GeoLib::GEOTYPE::POINT:
+		return this->getBoundaryElementsAtPoint(*dynamic_cast<const GeoLib::Point*>(&geoObj));
+		break;
 	case GeoLib::GEOTYPE::POLYLINE:
 		return this->getBoundaryElementsAlongPolyline(*dynamic_cast<const GeoLib::Polyline*>(&geoObj));
 		break;
@@ -49,6 +55,22 @@ std::vector<MeshLib::Element*> const& BoundaryElementsSearcher::getBoundaryEleme
 	}
 }
 
+std::vector<MeshLib::Element*> const&
+BoundaryElementsSearcher::getBoundaryElementsAtPoint(GeoLib::Point const& point)
+{
+	// look for already saved points and return if found.
+	for (auto const& boundaryElements : _boundary_elements_at_point)
+	{
+		if (boundaryElements->getPoint() == point)
+			return boundaryElements->getBoundaryElements();
+	}
+
+	// create new boundary elements at points.
+	_boundary_elements_at_point.push_back(
+	    new BoundaryElementsAtPoint(_mesh, _mshNodeSearcher, point));
+	return _boundary_elements_at_point.back()->getBoundaryElements();
+}
+
 std::vector<MeshLib::Element*> const& BoundaryElementsSearcher::getBoundaryElementsAlongPolyline(GeoLib::Polyline const& ply)
 {
 	std::vector<BoundaryElementsAlongPolyline*>::const_iterator it(_boundary_elements_along_polylines.begin());
diff --git a/MeshGeoToolsLib/BoundaryElementsSearcher.h b/MeshGeoToolsLib/BoundaryElementsSearcher.h
index 4e126f65dfa4db0192db6fac844123fe76440348..973bfee9045f54e0f80d93d6c71bf6d45bb822d0 100644
--- a/MeshGeoToolsLib/BoundaryElementsSearcher.h
+++ b/MeshGeoToolsLib/BoundaryElementsSearcher.h
@@ -13,6 +13,7 @@
 namespace GeoLib
 {
 class GeoObject;
+class Point;
 class Polyline;
 class Surface;
 }
@@ -26,6 +27,7 @@ class Element;
 namespace MeshGeoToolsLib
 {
 class MeshNodeSearcher;
+class BoundaryElementsAtPoint;
 class BoundaryElementsAlongPolyline;
 class BoundaryElementsOnSurface;
 
@@ -54,6 +56,13 @@ public:
 	 */
 	std::vector<MeshLib::Element*> const& getBoundaryElements(GeoLib::GeoObject const& geoObj);
 
+	/**
+	 * generate boundary elements at the given point.
+	 * @param point Search the mesh for given point
+	 * @return a vector of boundary elements
+	 */
+	std::vector<MeshLib::Element*> const& getBoundaryElementsAtPoint(
+	    GeoLib::Point const& point);
 	/**
 	 * generate boundary elements on the given polyline.
 	 * @param ply the GeoLib::Polyline the nearest mesh nodes are searched for
@@ -72,6 +81,7 @@ public:
 private:
 	MeshLib::Mesh const& _mesh;
 	MeshNodeSearcher &_mshNodeSearcher;
+	std::vector<BoundaryElementsAtPoint*> _boundary_elements_at_point;
 	std::vector<BoundaryElementsAlongPolyline*> _boundary_elements_along_polylines;
 	std::vector<BoundaryElementsOnSurface*> _boundary_elements_along_surfaces;
 };
diff --git a/MeshLib/Elements/Elements.h b/MeshLib/Elements/Elements.h
index 7e699cd0e3d99fc155a81f00d01df62c0811a994..b8e8395a370b63d6bd4aacb723c53cfe93d2fc42 100644
--- a/MeshLib/Elements/Elements.h
+++ b/MeshLib/Elements/Elements.h
@@ -16,6 +16,7 @@
 #include "MeshLib/Elements/Element.h"
 #include "MeshLib/Elements/Line.h"
 #include "MeshLib/Elements/Hex.h"
+#include "MeshLib/Elements/Point.h"
 #include "MeshLib/Elements/Prism.h"
 #include "MeshLib/Elements/Pyramid.h"
 #include "MeshLib/Elements/Quad.h"
diff --git a/MeshLib/Elements/Point.h b/MeshLib/Elements/Point.h
new file mode 100644
index 0000000000000000000000000000000000000000..02d50734c224e0d215e3e731e84ccb8fa128dbc1
--- /dev/null
+++ b/MeshLib/Elements/Point.h
@@ -0,0 +1,23 @@
+/**
+ * \copyright
+ * Copyright (c) 2012-2015, 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 MESHLIB_POINT_H_
+#define MESHLIB_POINT_H_
+
+#include "TemplateElement.h"
+#include "PointRule1.h"
+
+extern template class MeshLib::TemplateElement<MeshLib::PointRule1>;
+
+namespace MeshLib
+{
+	using Point = TemplateElement<PointRule1>;
+}
+
+#endif  // MESHLIB_POINT_H_
diff --git a/MeshLib/Elements/PointRule1.cpp b/MeshLib/Elements/PointRule1.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..c5f2bb68edaa72c11bfb2bf6beeee130ae5002f4
--- /dev/null
+++ b/MeshLib/Elements/PointRule1.cpp
@@ -0,0 +1,47 @@
+/**
+ * \copyright
+ * Copyright (c) 2012-2015, OpenGeoSys Community (http://www.opengeosys.org)
+ *            Distributed under a Modified BSD License.
+ *              See accompanying file LICENSE.txt or
+ *              http://www.opengeosys.org/project/license
+ *
+ */
+
+#include "PointRule1.h"
+
+#include "MathLib/Point3d.h"
+#include "MeshLib/Node.h"
+
+namespace MeshLib {
+
+const unsigned PointRule1::edge_nodes[1][1] =
+{
+	{0}
+};
+
+double PointRule1::computeVolume(Node const* const* /*_nodes*/)
+{
+	return 0;
+}
+
+bool PointRule1::isPntInElement(Node const* const* _nodes, MathLib::Point3d const& pnt, double eps)
+{
+	double const dist = MathLib::sqrDist(*_nodes[0], pnt);
+	return (dist < eps);
+}
+
+unsigned PointRule1::identifyFace(Node const* const* _nodes, Node* nodes[1])
+{
+	if (nodes[0] == _nodes[0])
+		return 0;
+	return std::numeric_limits<unsigned>::max();
+}
+
+ElementErrorCode PointRule1::validate(const Element* e)
+{
+	ElementErrorCode error_code;
+	error_code[ElementErrorFlag::ZeroVolume] = e->hasZeroVolume();
+	return error_code;
+}
+
+} // end namespace MeshLib
diff --git a/MeshLib/Elements/PointRule1.h b/MeshLib/Elements/PointRule1.h
new file mode 100644
index 0000000000000000000000000000000000000000..80da548d108d6ef03101729417291c66698044bc
--- /dev/null
+++ b/MeshLib/Elements/PointRule1.h
@@ -0,0 +1,70 @@
+/**
+ * \copyright
+ * Copyright (c) 2012-2015, 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 MESHLIB_POINTRULE1_H_
+#define MESHLIB_POINTRULE1_H_
+
+#include "MeshLib/MeshEnums.h"
+#include "Element.h"
+#include "VertexRule.h"
+#include "EdgeReturn.h"
+
+namespace MeshLib
+{
+
+/// A 0d point element.
+class PointRule1 : public VertexRule
+{
+public:
+	/// Constant: The number of base nodes for this element
+	static const unsigned n_base_nodes = 1u;
+
+	/// Constant: The number of all nodes for this element
+	static const unsigned n_all_nodes = 1u;
+
+	/// Constant: The geometric type of the element
+	static const MeshElemType mesh_elem_type = MeshElemType::POINT;
+
+	/// Constant: The FEM type of the element
+	static const CellType cell_type = CellType::POINT1;
+
+	/// Constant: The number of neighbors
+	static const unsigned n_neighbors = 2;
+
+	/// Constant: Local node index table for edge
+	static const unsigned edge_nodes[1][1];
+
+	/// Edge rule
+	typedef NoEdgeReturn EdgeReturn;
+
+	/// Checks if a point is inside the element.
+	///
+	/// Specifically this function tests if the squared euclidean distance
+	/// between the points \c _nodes and \c pnt is less then \c eps.
+	///
+	/// \param pnt a 3D MathLib::Point3d object
+	/// \param eps tolerance for numerical algorithm used for computing the
+	/// property
+	/// \return true if the point is not outside the element, false otherwise
+	static bool isPntInElement(Node const* const* _nodes,
+	                           MathLib::Point3d const& pnt, double eps);
+
+	/// Tests if the element is geometrically valid.
+	static ElementErrorCode validate(const Element* e);
+
+	/// Returns the ID of a face given an array of nodes.
+	static unsigned identifyFace(Node const* const*, Node* nodes[1]);
+
+	/// Calculates the length of a line
+	static double computeVolume(Node const* const* _nodes);
+};
+}
+
+#endif  // MESHLIB_POINTRULE1_H_
+
diff --git a/MeshLib/Elements/TemplateElement-impl.h b/MeshLib/Elements/TemplateElement-impl.h
index 7c9f88b5d57ad076db08f29300cc4adf6c95038c..cba91a51bb326b423673efeec0c369ee5ed735e3 100644
--- a/MeshLib/Elements/TemplateElement-impl.h
+++ b/MeshLib/Elements/TemplateElement-impl.h
@@ -46,13 +46,35 @@ TemplateElement<ELEMENT_RULE>::TemplateElement(const TemplateElement &e)
 	this->_content = e.getContent();
 }
 
+
+namespace detail
+{
+
+template<unsigned N>
+bool isEdge(unsigned const (&edge_nodes)[N], unsigned idx1, unsigned idx2)
+{
+
+	if (edge_nodes[0]==idx1 && edge_nodes[1]==idx2) return true;
+	if (edge_nodes[1]==idx1 && edge_nodes[0]==idx2) return true;
+
+	return false;
+}
+
+inline bool
+isEdge(unsigned const (&/*edge_nodes*/)[1], unsigned /*idx1*/, unsigned /*idx2*/)
+{
+	return false;
+}
+
+} // namespace detail
+
+
 template <class ELEMENT_RULE>
 bool TemplateElement<ELEMENT_RULE>::isEdge(unsigned idx1, unsigned idx2) const
 {
 	for (unsigned i(0); i<getNEdges(); i++)
 	{
-		if (ELEMENT_RULE::edge_nodes[i][0]==idx1 && ELEMENT_RULE::edge_nodes[i][1]==idx2) return true;
-		if (ELEMENT_RULE::edge_nodes[i][1]==idx1 && ELEMENT_RULE::edge_nodes[i][0]==idx2) return true;
+		if (detail::isEdge(ELEMENT_RULE::edge_nodes[i], idx1, idx2)) return true;
 	}
 	return false;
 }
diff --git a/MeshLib/Elements/TemplateElement.cpp b/MeshLib/Elements/TemplateElement.cpp
index 0934918e068b02610287e8e907dea8119543751e..dc07468f3d0facf2caced53c4c698134eeb69480 100644
--- a/MeshLib/Elements/TemplateElement.cpp
+++ b/MeshLib/Elements/TemplateElement.cpp
@@ -10,6 +10,7 @@
 #include "TemplateElement.h"
 
 #include "MeshLib/Elements/Hex.h"
+#include "MeshLib/Elements/Point.h"
 #include "MeshLib/Elements/Line.h"
 #include "MeshLib/Elements/Prism.h"
 #include "MeshLib/Elements/Pyramid.h"
@@ -33,6 +34,7 @@ template class MeshLib::TemplateElement<MeshLib::HexRule20>;
 template class MeshLib::TemplateElement<MeshLib::HexRule8>;
 template class MeshLib::TemplateElement<MeshLib::LineRule2>;
 template class MeshLib::TemplateElement<MeshLib::LineRule3>;
+template class MeshLib::TemplateElement<MeshLib::PointRule1>;
 template class MeshLib::TemplateElement<MeshLib::PrismRule15>;
 template class MeshLib::TemplateElement<MeshLib::PrismRule6>;
 template class MeshLib::TemplateElement<MeshLib::PyramidRule13>;
diff --git a/MeshLib/Elements/VertexRule.h b/MeshLib/Elements/VertexRule.h
new file mode 100644
index 0000000000000000000000000000000000000000..79a899172a36d355b153aced7ddcd5ce3c772e85
--- /dev/null
+++ b/MeshLib/Elements/VertexRule.h
@@ -0,0 +1,51 @@
+/**
+ * \copyright
+ * Copyright (c) 2012-2015, 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 MESHLIB_VERTEX_RULE_H
+#define MESHLIB_VERTEX_RULE_H
+
+namespace MeshLib
+{
+
+class Element;
+
+class VertexRule
+{
+public:
+	/// Constant: Dimension of this mesh element
+	static const unsigned dimension = 0u;
+
+	/// Constant: The number of faces
+	static const unsigned n_faces = 0;
+
+	/// Constant: The number of edges
+	static const unsigned n_edges = 0;
+
+	/// Get the number of nodes for face i.
+	static unsigned getNFaceNodes(unsigned /*i*/)
+	{
+		return 0;
+	}
+
+	/// Returns the i-th face of the element.
+	static const Element* getFace(const Element* /*e*/, unsigned /*i*/)
+	{
+		return nullptr;
+	}
+
+	/// Checks if the node order of an element is correct by testing surface
+	/// normals.  For 0D elements this always returns true.
+	static bool testElementNodeOrder(const Element* /*e*/) { return true; }
+
+};
+
+}
+
+#endif // MESHLIB_VERTEX_RULE_H
+
diff --git a/MeshLib/MeshEnums.cpp b/MeshLib/MeshEnums.cpp
index 76320d3646c69161f6791d93e1e47a7f803a0dde..b045cb3262a6577f8b415cacdf40092e96c73015 100644
--- a/MeshLib/MeshEnums.cpp
+++ b/MeshLib/MeshEnums.cpp
@@ -18,6 +18,8 @@ namespace MeshLib {
 
 const std::string MeshElemType2String(const MeshElemType t)
 {
+	if (t == MeshElemType::POINT)
+		return "Point";
 	if (t == MeshElemType::LINE)
 		return "Line";
 	if (t == MeshElemType::QUAD)
@@ -37,6 +39,8 @@ const std::string MeshElemType2String(const MeshElemType t)
 
 const std::string MeshElemType2StringShort(const MeshElemType t)
 {
+	if (t == MeshElemType::POINT)
+		return "point";
 	if (t == MeshElemType::LINE)
 		return "line";
 	if (t == MeshElemType::QUAD)
@@ -56,6 +60,8 @@ const std::string MeshElemType2StringShort(const MeshElemType t)
 
 MeshElemType String2MeshElemType(const std::string &s)
 {
+	if ((s.compare("point") == 0) || (s.compare("Point") == 0))
+		return MeshElemType::POINT;
 	if ((s.compare("line") == 0) || (s.compare("Line") == 0))
 		return MeshElemType::LINE;
 	if ((s.compare("quad") == 0) || (s.compare("Quadrilateral") == 0))
@@ -76,6 +82,7 @@ MeshElemType String2MeshElemType(const std::string &s)
 std::vector<MeshElemType> getMeshElemTypes()
 {
 	std::vector<MeshElemType> vec;
+	vec.push_back(MeshElemType::POINT);
 	vec.push_back(MeshElemType::LINE);
 	vec.push_back(MeshElemType::QUAD);
 	vec.push_back(MeshElemType::HEXAHEDRON);
@@ -100,6 +107,7 @@ const std::string CellType2String(const CellType t)
 	if (t == CellType::type)\
 		return #type;
 
+	RETURN_CELL_TYPE_STR(t, POINT1);
 	RETURN_CELL_TYPE_STR(t, LINE2);
 	RETURN_CELL_TYPE_STR(t, LINE3);
 	RETURN_CELL_TYPE_STR(t, QUAD4);
diff --git a/MeshLib/MeshEnums.h b/MeshLib/MeshEnums.h
index 3a08a551d6b7a27f4e66083af60079f51db52fc1..dfa775ef8f49a2a535d17f34202ca9061bb3fce5 100644
--- a/MeshLib/MeshEnums.h
+++ b/MeshLib/MeshEnums.h
@@ -27,6 +27,7 @@ namespace MeshLib {
 enum class MeshElemType
 {
 	INVALID = 0,
+	POINT = 1,
 	LINE = 3,
 	QUAD = 9,
 	HEXAHEDRON = 12,
@@ -42,6 +43,7 @@ enum class MeshElemType
 enum class CellType
 {
 	INVALID,
+	POINT1,
 	LINE2,
 	LINE3,
 	TRI3,
diff --git a/NumLib/Fem/CoordinatesMapping/NaturalCoordinatesMapping-impl.h b/NumLib/Fem/CoordinatesMapping/NaturalCoordinatesMapping-impl.h
index f657cd7259b4b2d0108d5e04c1304fa2a652f994..207bc7bec0663c4c0fad8fd3de229d793ed8d7f5 100644
--- a/NumLib/Fem/CoordinatesMapping/NaturalCoordinatesMapping-impl.h
+++ b/NumLib/Fem/CoordinatesMapping/NaturalCoordinatesMapping-impl.h
@@ -36,7 +36,9 @@ inline void computeMappingMatrices(
 }
 
 template <class T_MESH_ELEMENT, class T_SHAPE_FUNC, class T_SHAPE_MATRICES>
-inline void computeMappingMatrices(
+inline
+typename std::enable_if<T_SHAPE_FUNC::DIM!=0>::type
+computeMappingMatrices(
         const T_MESH_ELEMENT &/*ele*/,
         const double* natural_pt,
         const MeshLib::ElementCoordinatesMappingLocal &/*ele_local_coord*/,
@@ -46,9 +48,22 @@ inline void computeMappingMatrices(
     double* const dNdr = shapemat.dNdr.data();
     T_SHAPE_FUNC::computeGradShapeFunction(natural_pt, dNdr);
 }
+template <class T_MESH_ELEMENT, class T_SHAPE_FUNC, class T_SHAPE_MATRICES>
+inline
+typename std::enable_if<T_SHAPE_FUNC::DIM==0>::type
+computeMappingMatrices(
+        const T_MESH_ELEMENT &/*ele*/,
+        const double* /*natural_pt*/,
+        const MeshLib::ElementCoordinatesMappingLocal &/*ele_local_coord*/,
+        T_SHAPE_MATRICES &/*shapemat*/,
+        FieldType<ShapeMatrixType::DNDR>)
+{
+}
 
 template <class T_MESH_ELEMENT, class T_SHAPE_FUNC, class T_SHAPE_MATRICES>
-inline void computeMappingMatrices(
+inline
+typename std::enable_if<T_SHAPE_FUNC::DIM!=0>::type
+computeMappingMatrices(
         const T_MESH_ELEMENT &ele,
         const double* natural_pt,
         const MeshLib::ElementCoordinatesMappingLocal &ele_local_coord,
@@ -79,6 +94,18 @@ inline void computeMappingMatrices(
         ERR("***error: det|J|=%e is not positive.\n", shapemat.detJ);
 #endif
 }
+template <class T_MESH_ELEMENT, class T_SHAPE_FUNC, class T_SHAPE_MATRICES>
+inline
+typename std::enable_if<T_SHAPE_FUNC::DIM==0>::type
+computeMappingMatrices(
+        const T_MESH_ELEMENT &/*ele*/,
+        const double* /*natural_pt*/,
+        const MeshLib::ElementCoordinatesMappingLocal &/*ele_local_coord*/,
+        T_SHAPE_MATRICES &shapemat,
+        FieldType<ShapeMatrixType::DNDR_J>)
+{
+    shapemat.detJ = 1.0;
+}
 
 template <class T_MESH_ELEMENT, class T_SHAPE_FUNC, class T_SHAPE_MATRICES>
 inline void computeMappingMatrices(
@@ -95,7 +122,9 @@ inline void computeMappingMatrices(
 }
 
 template <class T_MESH_ELEMENT, class T_SHAPE_FUNC, class T_SHAPE_MATRICES>
-inline void computeMappingMatrices(
+inline
+typename std::enable_if<T_SHAPE_FUNC::DIM!=0>::type
+computeMappingMatrices(
         const T_MESH_ELEMENT &ele,
         const double* natural_pt,
         const MeshLib::ElementCoordinatesMappingLocal &ele_local_coord,
@@ -124,6 +153,20 @@ inline void computeMappingMatrices(
         }
     }
 }
+template <class T_MESH_ELEMENT, class T_SHAPE_FUNC, class T_SHAPE_MATRICES>
+inline
+typename std::enable_if<T_SHAPE_FUNC::DIM==0>::type
+computeMappingMatrices(
+       const T_MESH_ELEMENT &ele,
+       const double* natural_pt,
+       const MeshLib::ElementCoordinatesMappingLocal &ele_local_coord,
+       T_SHAPE_MATRICES &shapemat,
+       FieldType<ShapeMatrixType::DNDX>)
+{
+    computeMappingMatrices<T_MESH_ELEMENT, T_SHAPE_FUNC, T_SHAPE_MATRICES>
+        (ele, natural_pt, ele_local_coord, shapemat, FieldType<ShapeMatrixType::DNDR_J>());
+}
+
 
 template <class T_MESH_ELEMENT, class T_SHAPE_FUNC, class T_SHAPE_MATRICES>
 inline void computeMappingMatrices(
diff --git a/NumLib/Fem/FiniteElement/C0IsoparametricElements.h b/NumLib/Fem/FiniteElement/C0IsoparametricElements.h
index 7da4ddc6ac1477837bc1a0f38c37ec2d0d85023a..3cc66d3a5c5db39833b783c1247f8c58426635da 100644
--- a/NumLib/Fem/FiniteElement/C0IsoparametricElements.h
+++ b/NumLib/Fem/FiniteElement/C0IsoparametricElements.h
@@ -14,6 +14,7 @@
 #ifndef C0ISOPARAMETRICELEMENTS_H_
 #define C0ISOPARAMETRICELEMENTS_H_
 
+#include "NumLib/Fem/ShapeFunction/ShapePoint1.h"
 #include "NumLib/Fem/ShapeFunction/ShapeLine2.h"
 #include "NumLib/Fem/ShapeFunction/ShapeLine3.h"
 #include "NumLib/Fem/ShapeFunction/ShapeTri3.h"
@@ -35,6 +36,12 @@
 namespace NumLib
 {
 
+template <template <typename> class T_SHAPE_MATRIX_POLICY>
+struct FePOINT1
+{
+    typedef TemplateIsoparametric<ShapePoint1, T_SHAPE_MATRIX_POLICY<ShapePoint1>> type;
+};
+
 template <template <typename> class T_SHAPE_MATRIX_POLICY>
 struct FeLINE2
 {
diff --git a/NumLib/Fem/Integration/GaussIntegrationPolicy.h b/NumLib/Fem/Integration/GaussIntegrationPolicy.h
index 937a38ac62c96b9bb0606bf3db8b2d3702339428..5503e0f13c7658f5ab42dddbd4fbdbc3dba127f6 100644
--- a/NumLib/Fem/Integration/GaussIntegrationPolicy.h
+++ b/NumLib/Fem/Integration/GaussIntegrationPolicy.h
@@ -10,6 +10,7 @@
 #ifndef NUMLIB_GAUSSINTEGRATIONPOLICY_H_
 #define NUMLIB_GAUSSINTEGRATIONPOLICY_H_
 
+#include "MeshLib/Elements/Point.h"
 #include "MeshLib/Elements/Tri.h"
 #include "MeshLib/Elements/Tet.h"
 #include "MeshLib/Elements/Prism.h"
@@ -20,6 +21,7 @@
 #include "NumLib/Fem/Integration/IntegrationGaussTet.h"
 #include "NumLib/Fem/Integration/IntegrationGaussPrism.h"
 #include "NumLib/Fem/Integration/IntegrationGaussPyramid.h"
+#include "NumLib/Fem/Integration/IntegrationPoint.h"
 
 namespace NumLib
 {
@@ -40,6 +42,13 @@ struct GaussIntegrationPolicy
         NumLib::IntegrationGaussRegular<MeshElement::dimension>;
 };
 
+template<>
+struct GaussIntegrationPolicy<MeshLib::Point>
+{
+    using MeshElement = MeshLib::Point;
+    using IntegrationMethod = NumLib::IntegrationPoint;
+};
+
 template<>
 struct GaussIntegrationPolicy<MeshLib::Tri>
 {
diff --git a/NumLib/Fem/Integration/IntegrationPoint.h b/NumLib/Fem/Integration/IntegrationPoint.h
new file mode 100644
index 0000000000000000000000000000000000000000..056b69ebf26ba4eab4e88cc5d83736f347d3e296
--- /dev/null
+++ b/NumLib/Fem/Integration/IntegrationPoint.h
@@ -0,0 +1,83 @@
+/**
+ * \copyright
+ * Copyright (c) 2012-2015, 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 NUMLIB_INTEGRATIONPOINT_H_
+#define NUMLIB_INTEGRATIONPOINT_H_
+
+namespace NumLib
+{
+/// Integration rule for point elements.
+///
+/// The integration order is not stored or used for point integration.
+/// It is only needed to satisfy the common integration rule concepts.
+class IntegrationPoint
+{
+	typedef MathLib::TemplateWeightedPoint<double, double, 1> WeightedPoint;
+
+public:
+	/// IntegrationPoint constructor for given order.
+	explicit IntegrationPoint(std::size_t /* order */)
+	{
+	}
+
+	/// Change the integration order.
+	void setIntegrationOrder(std::size_t /* order */)
+	{
+	}
+
+	/// Return current integration order.
+	std::size_t getIntegrationOrder() const
+	{
+		return 0;
+	}
+
+	/// Return the number of sampling points.
+	std::size_t getNPoints() const
+	{
+		return 1;
+	}
+
+	/// Get coordinates of a integration point.
+	///
+	/// \param igp      The integration point index.
+	/// \return a weighted point.
+	WeightedPoint getWeightedPoint(std::size_t igp)
+	{
+		return getWeightedPoint(getIntegrationOrder(), igp);
+	}
+
+	/// Get coordinates of a integration point.
+	///
+	/// \param order    the number of integration points.
+	/// \param pt_id     the sampling point id.
+	/// \return weight
+	static WeightedPoint getWeightedPoint(std::size_t /* order */,
+	                                      std::size_t /*igp*/)
+	{
+		return WeightedPoint({{1}}, 1);
+	}
+
+	template <typename Method>
+	static WeightedPoint getWeightedPoint(std::size_t /*igp*/)
+	{
+		return WeightedPoint({{1}}, 1);
+	}
+
+	/// Get the number of integration points.
+	///
+	/// \param order    the number of integration points
+	/// \return the number of points.
+	static std::size_t getNPoints(std::size_t /* order */)
+	{
+		return 1;
+	}
+};
+}
+
+#endif  // NUMLIB_INTEGRATIONPOINT_H_
diff --git a/NumLib/Fem/ShapeFunction/ShapePoint1-impl.h b/NumLib/Fem/ShapeFunction/ShapePoint1-impl.h
new file mode 100644
index 0000000000000000000000000000000000000000..81907f2162d79b359fed5f6bb2e87195f656b31e
--- /dev/null
+++ b/NumLib/Fem/ShapeFunction/ShapePoint1-impl.h
@@ -0,0 +1,19 @@
+/**
+ * \copyright
+ * Copyright (c) 2012-2015, 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 ShapePoint1::computeShapeFunction(const T_X& /*r*/, T_N& N)
+{
+	N[0] = 1;
+}
+
+}
diff --git a/NumLib/Fem/ShapeFunction/ShapePoint1.h b/NumLib/Fem/ShapeFunction/ShapePoint1.h
new file mode 100644
index 0000000000000000000000000000000000000000..8ce677719c86e0347d0c0803b764d1b0d5e99790
--- /dev/null
+++ b/NumLib/Fem/ShapeFunction/ShapePoint1.h
@@ -0,0 +1,36 @@
+/**
+ * \copyright
+ * Copyright (c) 2012-2015, 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 SHAPEPOINT1_H_
+#define SHAPEPOINT1_H_
+
+#include "MeshLib/Elements/Point.h"
+
+namespace NumLib
+{
+///  Shape function for a point element in natural coordinates.
+class ShapePoint1
+{
+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);
+
+	using MeshElement = MeshLib::Point;
+	static const std::size_t DIM = MeshElement::dimension;
+	static const std::size_t NPOINTS = MeshElement::n_all_nodes;
+};
+}
+
+#include "ShapePoint1-impl.h"
+
+#endif  // SHAPEPOINT1_H_
diff --git a/NumLib/Fem/ShapeMatrixPolicy.h b/NumLib/Fem/ShapeMatrixPolicy.h
index 5e8c9ae73b6a467756e5f814386036a27c4169c9..2146bba64703774b26fe1ce4201936d18814313a 100644
--- a/NumLib/Fem/ShapeMatrixPolicy.h
+++ b/NumLib/Fem/ShapeMatrixPolicy.h
@@ -13,17 +13,54 @@
 #include "NumLib/Fem/CoordinatesMapping/ShapeMatrices.h"
 
 #ifdef OGS_USE_EIGEN
-#include <Eigen/Eigen>
+#include <Eigen/Dense>
+
+namespace detail
+{
+    /// Forwards the Eigen::Matrix type for general N and M.
+    /// There is a partial specialization for M = 1 to store the matrix in
+    /// column major storage order.
+    template <int N, int M>
+    struct EigenMatrixType
+    {
+        using type = Eigen::Matrix<double, N, M, Eigen::RowMajor>;
+    };
+
+    /// Specialization for Nx1 matrices which can be stored only in column major
+    /// form in Eigen-3.2.5.
+    template <int N>
+    struct EigenMatrixType<N, 1>
+    {
+        using type = Eigen::Matrix<double, N, 1, Eigen::ColMajor>;
+    };
+
+    /// Specialization for 0xM matrices. Using fixed size Eigen matrices here
+    /// would lead to zero sized arrays, which cause compilation errors on
+    /// some compilers.
+    template <int M>
+    struct EigenMatrixType<0, M>
+    {
+        using type = Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>;
+    };
+
+    template <>
+    struct EigenMatrixType<0, 1>
+    {
+        using type = Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic, Eigen::ColMajor>;
+    };
+
+}   // detail
 
 /// An implementation of ShapeMatrixPolicy using fixed size (compile-time) eigen
 /// matrices and vectors.
 template <typename ShapeFunction, unsigned GlobalDim>
 struct EigenFixedShapeMatrixPolicy
 {
-    template <int N, int M>
-    using _MatrixType = Eigen::Matrix<double, N, M, Eigen::RowMajor>;
     template <int N>
-    using _VectorType = Eigen::Matrix<double, N, 1>;
+    using _VectorType = typename ::detail::EigenMatrixType<N, 1>::type;
+
+    template <int N, int M>
+    using _MatrixType = typename ::detail::EigenMatrixType<N, M>::type;
 
     using NodalMatrixType = _MatrixType<ShapeFunction::NPOINTS, ShapeFunction::NPOINTS>;
     using NodalVectorType = _VectorType<ShapeFunction::NPOINTS>;
diff --git a/Tests/Data/Elliptic/line_1_GroundWaterFlow/line_1.gml.md5 b/Tests/Data/Elliptic/line_1_GroundWaterFlow/line_1.gml.md5
new file mode 100644
index 0000000000000000000000000000000000000000..2fdd2fda0d910aee6737f0e716c12c92f6ffe0c2
--- /dev/null
+++ b/Tests/Data/Elliptic/line_1_GroundWaterFlow/line_1.gml.md5
@@ -0,0 +1 @@
+4ec0050b42a34e66c9fd5eeca418e1e0
diff --git a/Tests/Data/Elliptic/line_1_GroundWaterFlow/line_1_line_1e1.vtu.md5 b/Tests/Data/Elliptic/line_1_GroundWaterFlow/line_1_line_1e1.vtu.md5
new file mode 100644
index 0000000000000000000000000000000000000000..e0b70e052cc1961e5e63933c6b6d7c05f3405210
--- /dev/null
+++ b/Tests/Data/Elliptic/line_1_GroundWaterFlow/line_1_line_1e1.vtu.md5
@@ -0,0 +1 @@
+043336d543300ea41fd5afd8412df4db
diff --git a/Tests/Data/Elliptic/line_1_GroundWaterFlow/line_1e1.prj.md5 b/Tests/Data/Elliptic/line_1_GroundWaterFlow/line_1e1.prj.md5
new file mode 100644
index 0000000000000000000000000000000000000000..2b639a55d9d41a769a6ade20c7609c675d2340c7
--- /dev/null
+++ b/Tests/Data/Elliptic/line_1_GroundWaterFlow/line_1e1.prj.md5
@@ -0,0 +1 @@
+250962ebd5a5b4fb0ef91fd3ab19b465
diff --git a/Tests/Data/Elliptic/line_1_GroundWaterFlow/line_1e1_neumann.prj.md5 b/Tests/Data/Elliptic/line_1_GroundWaterFlow/line_1e1_neumann.prj.md5
new file mode 100644
index 0000000000000000000000000000000000000000..0fe415be3c42e709c4392d4d09c36b3016c6f7f7
--- /dev/null
+++ b/Tests/Data/Elliptic/line_1_GroundWaterFlow/line_1e1_neumann.prj.md5
@@ -0,0 +1 @@
+7dab35ae99f9d34d53ebe42d5ba7afaf