added angle-parameter to extract surfaces more precisely

Applications/ApplicationsLib/ProjectData.cpp

@@ -84,12 +84,14 @@ void ProjectData::addMesh(MeshLib::Mesh* mesh)
 	_mesh_vec.push_back(mesh);
 }
 
+
 std::vector<MeshLib::Mesh*>::const_iterator ProjectData::findMeshByName(
 		std::string const& name) const
 {
 	return const_cast<ProjectData&>(*this).findMeshByName(name);
 }
 
+
 std::vector<MeshLib::Mesh*>::iterator ProjectData::findMeshByName(
 		std::string const& name)
 {
@@ -98,7 +100,6 @@ std::vector<MeshLib::Mesh*>::iterator ProjectData::findMeshByName(
 			{
 				return mesh && (name == mesh->getName());
 			});
-
 }
 
 const MeshLib::Mesh* ProjectData::getMesh(const std::string &name) const

Applications/DataExplorer/DataView/DirectConditionGenerator.cpp

@@ -41,7 +41,7 @@ const std::vector< std::pair<size_t,double> >& DirectConditionGenerator::directT
 		}
 
 		const MathLib::Vector3 dir(0,0,-1);
-		const std::vector<GeoLib::PointWithID*> surface_nodes(MeshLib::MeshSurfaceExtraction::getSurfaceNodes(mesh, dir) );
+		const std::vector<GeoLib::PointWithID*> surface_nodes(MeshLib::MeshSurfaceExtraction::getSurfaceNodes(mesh, dir, 90) );
 		const size_t nNodes(surface_nodes.size());
 		const double no_data (raster->getNoDataValue());
 		_direct_values.reserve(nNodes);

Applications/DataExplorer/DataView/GeoMapper.cpp

@@ -68,7 +68,7 @@ void GeoMapper::mapOnMesh(const MeshLib::Mesh* mesh)
 	else
 	{
 		const MathLib::Vector3 dir(0,0,-1);
-		this->_mesh = MeshLib::MeshSurfaceExtraction::getMeshSurface(*mesh, dir);
+		this->_mesh = MeshLib::MeshSurfaceExtraction::getMeshSurface(*mesh, dir, 90);
 	}
 	std::vector<GeoLib::PointWithID*> sfc_pnts;
 	// init grid

Applications/DataExplorer/DataView/MshView.cpp

@@ -194,7 +194,7 @@ void MshView::extractSurfaceMesh()
 
 	const MeshLib::Mesh* mesh = static_cast<MshModel*>(this->model())->getMesh(index);
 	const MathLib::Vector3 dir(0, 0, -1);
-	static_cast<MshModel*>(this->model())->addMesh( MeshLib::MeshSurfaceExtraction::getMeshSurface(*mesh, dir) );
+	static_cast<MshModel*>(this->model())->addMesh( MeshLib::MeshSurfaceExtraction::getMeshSurface(*mesh, dir, 90) );
 }
 
 void MshView::convertMeshToGeometry()

MathLib/Vector3.h

@@ -115,6 +115,14 @@ public:
 			this->_x[i] *= s;
 	}
 
+	/// Returns a normalized version of this vector
+	TemplateVector3<double> normalizeVector() const
+	{
+		TemplateVector3<double> norm_vec (_x[0], _x[1], _x[2]);
+		norm_vec.normalize();
+		return norm_vec;
+	}
+
 	/// Returns the squared length
 	double getSqrLength(void) const
 	{

MeshLib/Elements/Cell.cpp

@@ -36,11 +36,7 @@ Cell::~Cell()
 
 bool Cell::isOnSurface() const
 {
-	unsigned n (this->getNNeighbors());
-	for (unsigned i(0); i<n; i++)
-		if (!this->_neighbors[i])
-			return true;
-	return false;
+	return isBoundaryElement();
 }
 
 bool Cell::testElementNodeOrder() const

MeshLib/MeshSurfaceExtraction.cpp

@@ -15,6 +15,7 @@
 #include "MeshSurfaceExtraction.h"
 
 #include <cassert>
+#include <math.h>
 
 #include "logog/include/logog.hpp"
 
@@ -63,14 +64,17 @@ void MeshSurfaceExtraction::getSurfaceAreaForNodes(const MeshLib::Mesh &mesh, st
 		ERR ("Error in MeshSurfaceExtraction::getSurfaceAreaForNodes() - Given mesh is no surface mesh (dimension != 2).");
 }
 
-MeshLib::Mesh* MeshSurfaceExtraction::getMeshSurface(const MeshLib::Mesh &mesh, const MathLib::Vector3 &dir, bool keepOriginalNodeIds)
+MeshLib::Mesh* MeshSurfaceExtraction::getMeshSurface(const MeshLib::Mesh &mesh, const MathLib::Vector3 &dir, double angle, bool keepOriginalNodeIds)
 {
+	if (angle< 0) angle = 0;
+	if (angle > 90) angle = 90;
+
 	INFO ("Extracting mesh surface...");
 	const std::vector<MeshLib::Element*> all_elements (mesh.getElements());
 	const std::vector<MeshLib::Node*> all_nodes (mesh.getNodes());
 
 	std::vector<MeshLib::Element*> sfc_elements;
-	get2DSurfaceElements(all_elements, sfc_elements, dir, mesh.getDimension());
+	get2DSurfaceElements(all_elements, sfc_elements, dir, angle, mesh.getDimension());
 
 	if (sfc_elements.empty())
 		return nullptr;
@@ -114,7 +118,7 @@ MeshLib::Mesh* MeshSurfaceExtraction::getMeshSurface(const MeshLib::Mesh &mesh,
 	return result;
 }
 
-void MeshSurfaceExtraction::get2DSurfaceElements(const std::vector<MeshLib::Element*> &all_elements, std::vector<MeshLib::Element*> &sfc_elements, const MathLib::Vector3 &dir, unsigned mesh_dimension)
+void MeshSurfaceExtraction::get2DSurfaceElements(const std::vector<MeshLib::Element*> &all_elements, std::vector<MeshLib::Element*> &sfc_elements, const MathLib::Vector3 &dir, double angle, unsigned mesh_dimension)
 {
 	if (mesh_dimension<2 || mesh_dimension>3)
 		ERR("Cannot handle meshes of dimension %i", mesh_dimension);
@@ -123,6 +127,10 @@ void MeshSurfaceExtraction::get2DSurfaceElements(const std::vector<MeshLib::Elem
 	if (MathLib::scalarProduct(dir, dir) != 0)
 		complete_surface = false;
 
+	double const pi (3.14159265358979323846);
+	double const cos_theta (std::cos(angle * pi / 180.0));
+	MathLib::Vector3 const norm_dir (dir.normalizeVector());
+
 	for (auto elem = all_elements.begin(); elem != all_elements.end(); ++elem)
 	{
 		const unsigned element_dimension ((*elem)->getDimension());
@@ -134,8 +142,7 @@ void MeshSurfaceExtraction::get2DSurfaceElements(const std::vector<MeshLib::Elem
 			if (!complete_surface)
 			{
 				MeshLib::Face* face = dynamic_cast<MeshLib::Face*>(*elem);
-
-				if (MathLib::scalarProduct(face->getSurfaceNormal(), dir) <= 0)
+				if (MathLib::scalarProduct(face->getSurfaceNormal().normalizeVector(), norm_dir) >= cos_theta)
 					continue;	
 			}
 			sfc_elements.push_back(*elem);
@@ -153,11 +160,13 @@ void MeshSurfaceExtraction::get2DSurfaceElements(const std::vector<MeshLib::Elem
 
 				const MeshLib::Face* face = static_cast<const MeshLib::Face*>(cell->getFace(j));
 				if (!complete_surface)
-					if (MathLib::scalarProduct(face->getSurfaceNormal(), dir) <= 0)
+				{
+					if (MathLib::scalarProduct(face->getSurfaceNormal().normalizeVector(), norm_dir) < cos_theta)
 					{
 						delete face;
 						continue;
 					}
+				}
 				if (face->getGeomType() == MeshElemType::TRIANGLE)
 					sfc_elements.push_back(new MeshLib::Tri(*static_cast<const MeshLib::Tri*>(face)));
 				else
@@ -191,14 +200,14 @@ void MeshSurfaceExtraction::get2DSurfaceNodes(const std::vector<MeshLib::Node*>
 	}
 }
 
-std::vector<GeoLib::PointWithID*> MeshSurfaceExtraction::getSurfaceNodes(const MeshLib::Mesh &mesh, const MathLib::Vector3 &dir)
+std::vector<GeoLib::PointWithID*> MeshSurfaceExtraction::getSurfaceNodes(const MeshLib::Mesh &mesh, const MathLib::Vector3 &dir, double angle)
 {
 	INFO ("Extracting surface nodes...");
 	const std::vector<MeshLib::Element*> all_elements (mesh.getElements());
 	const std::vector<MeshLib::Node*> all_nodes (mesh.getNodes());
 
 	std::vector<MeshLib::Element*> sfc_elements;
-	get2DSurfaceElements(all_elements, sfc_elements, dir, mesh.getDimension());
+	get2DSurfaceElements(all_elements, sfc_elements, dir, angle, mesh.getDimension());
 
 	std::vector<MeshLib::Node*> sfc_nodes;
 	std::vector<unsigned> node_id_map(mesh.getNNodes());

MeshLib/MeshSurfaceExtraction.h

@@ -40,20 +40,21 @@ public:
 	static void getSurfaceAreaForNodes(const MeshLib::Mesh &mesh, std::vector<double> &node_area_vec);
 
 	/// Returns the surface nodes of a layered mesh.
-	static std::vector<GeoLib::PointWithID*> getSurfaceNodes(const MeshLib::Mesh &mesh, const MathLib::Vector3 &dir);
+	static std::vector<GeoLib::PointWithID*> getSurfaceNodes(const MeshLib::Mesh &mesh, const MathLib::Vector3 &dir, double angle);
 
 	/**
 	 * Returns the 2d-element mesh representing the surface of the given layered mesh.
 	 * \param mesh                The original mesh
 	 * \param dir                 The direction in which face normals have to point to be considered surface elements
+     * \param angle               The angle of the allowed deviation from the given direction (0 < angle < 90 degrees)
 	 * \param keepOriginalNodeIds If true, ids of mesh nodes are set to ids in original mesh, otherwise node ids are reset (as usual when creating a mesh)
 	 * \return                    A 2D mesh representing the surface in direction dir
 	 */
-	static MeshLib::Mesh* getMeshSurface(const MeshLib::Mesh &mesh, const MathLib::Vector3 &dir, bool keepOriginalNodeIds = false);
+	static MeshLib::Mesh* getMeshSurface(const MeshLib::Mesh &mesh, const MathLib::Vector3 &dir, double angle, bool keepOriginalNodeIds = false);
 
 private:
 	/// Functionality needed for getSurfaceNodes() and getMeshSurface()
-	static void get2DSurfaceElements(const std::vector<MeshLib::Element*> &all_elements, std::vector<MeshLib::Element*> &sfc_elements, const MathLib::Vector3 &dir, unsigned mesh_dimension);
+	static void get2DSurfaceElements(const std::vector<MeshLib::Element*> &all_elements, std::vector<MeshLib::Element*> &sfc_elements, const MathLib::Vector3 &dir, double angle, unsigned mesh_dimension);
 
 	/// Functionality needed for getSurfaceNodes() and getMeshSurface()
 	static void get2DSurfaceNodes(const std::vector<MeshLib::Node*> &all_nodes, std::vector<MeshLib::Node*> &sfc_nodes, const std::vector<MeshLib::Element*> &sfc_elements, std::vector<unsigned> &node_id_map);

MeshLib/Node.h

@@ -37,7 +37,7 @@ class Node : public GeoLib::PointWithID
 {
 	/* friend functions: */
 	friend bool MeshLayerMapper::layerMapping(MeshLib::Mesh &mesh, const GeoLib::Raster &raster, double noDataReplacementValue);
-	friend MeshLib::Mesh* MeshSurfaceExtraction::getMeshSurface(const MeshLib::Mesh &mesh, const MathLib::Vector3 &dir, bool keep3dMeshIds);
+	friend MeshLib::Mesh* MeshSurfaceExtraction::getMeshSurface(const MeshLib::Mesh &mesh, const MathLib::Vector3 &dir, double angle, bool keep3dMeshIds);
 
 	/* friend classes: */
 	friend class Mesh;