Merge pull request #873 from chleh/dima-0dElement

0d element

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()

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))] =

MeshGeoToolsLib/BoundaryElementsAtPoint.cpp

+/**
+ * @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

MeshGeoToolsLib/BoundaryElementsAtPoint.h

+/**
+ * @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_

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());

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;
 };

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"

MeshLib/Elements/Point.h

+/**
+ * \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_

MeshLib/Elements/PointRule1.cpp

+/**
+ * \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

MeshLib/Elements/PointRule1.h

+/**
+ * \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_
+

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;
 }

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>;

MeshLib/Elements/VertexRule.h

+/**
+ * \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
+

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);

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,

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(

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
 {

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>
 {

NumLib/Fem/Integration/IntegrationPoint.h

+/**
+ * \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_

NumLib/Fem/ShapeFunction/ShapePoint1-impl.h

+/**
+ * \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;
+}
+
+}