diff --git a/Applications/ApplicationsLib/ProjectData.cpp b/Applications/ApplicationsLib/ProjectData.cpp index d0d5e2f9fbc007d086e5670d370ba9a84c5d7922..5ba3c4933e6860e449c07c870e23741d5fae4b91 100644 --- a/Applications/ApplicationsLib/ProjectData.cpp +++ b/Applications/ApplicationsLib/ProjectData.cpp @@ -98,7 +98,6 @@ std::vector<MeshLib::Mesh*>::iterator ProjectData::findMeshByName( { return mesh && (name == mesh->getName()); }); - } const MeshLib::Mesh* ProjectData::getMesh(const std::string &name) const diff --git a/Applications/DataExplorer/DataView/DirectConditionGenerator.cpp b/Applications/DataExplorer/DataView/DirectConditionGenerator.cpp index 89d63df52262a81c0c7cae3fdc7735c86e1baa1e..ec002fb6e7e14c415540792b9cc7d9d59a38a3ec 100644 --- a/Applications/DataExplorer/DataView/DirectConditionGenerator.cpp +++ b/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); @@ -72,8 +72,8 @@ const std::vector< std::pair<size_t,double> >& DirectConditionGenerator::directW const MathLib::Vector3 dir(0,0,-1); MeshLib::Mesh* sfc_mesh (MeshLib::MeshSurfaceExtraction::getMeshSurface(mesh, dir, true)); - std::vector<double> node_area_vec; - MeshLib::MeshSurfaceExtraction::getSurfaceAreaForNodes(*sfc_mesh, node_area_vec); + std::vector<double> node_area_vec = + MeshLib::MeshSurfaceExtraction::getSurfaceAreaForNodes(*sfc_mesh); const std::vector<MeshLib::Node*> &surface_nodes (sfc_mesh->getNodes()); const size_t nNodes(sfc_mesh->getNNodes()); const double no_data (raster->getNoDataValue()); diff --git a/Applications/DataExplorer/DataView/GeoMapper.cpp b/Applications/DataExplorer/DataView/GeoMapper.cpp index 9541cb730bdcb41fe430bcc96d0507135253e0f8..9684ea30e3bb86decf4e447d3896b22feb9aef13 100644 --- a/Applications/DataExplorer/DataView/GeoMapper.cpp +++ b/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 diff --git a/Applications/DataExplorer/DataView/MshView.cpp b/Applications/DataExplorer/DataView/MshView.cpp index c7f512a66ebc966acf804c093352b916d18e20f2..79cda6174da0f17e31721f2ad427f0b1fdacf6e0 100644 --- a/Applications/DataExplorer/DataView/MshView.cpp +++ b/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, 89) ); } void MshView::convertMeshToGeometry() diff --git a/MathLib/Vector3.h b/MathLib/Vector3.h index d218249698a0cb72c4384a3dbc8e384de5cc47cf..5fc54e6918289ed0e35ee3b2bafcf78c6114107a 100644 --- a/MathLib/Vector3.h +++ b/MathLib/Vector3.h @@ -115,6 +115,14 @@ public: this->_x[i] *= s; } + /// Returns a normalized version of this vector + TemplateVector3<double> getNormalizedVector() const + { + TemplateVector3<double> norm_vec (this->_x[0], this->_x[1], this->_x[2]); + norm_vec.normalize(); + return norm_vec; + } + /// Returns the squared length double getSqrLength(void) const { diff --git a/MeshLib/Elements/Cell.cpp b/MeshLib/Elements/Cell.cpp index 4bb9c7f9a13f686ba5446cbae338502eb9f9850b..45d5dca6e9be13342aec9957ff1a608baa9c1855 100644 --- a/MeshLib/Elements/Cell.cpp +++ b/MeshLib/Elements/Cell.cpp @@ -34,15 +34,6 @@ Cell::Cell(unsigned value, std::size_t id) 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; -} - bool Cell::testElementNodeOrder() const { const MathLib::Vector3 c (getCenterOfGravity()); diff --git a/MeshLib/Elements/Cell.h b/MeshLib/Elements/Cell.h index 62b1410abcc654fdc8ead5110e59a1d3d4ca7db1..dfe4b9b9946df75c0f6c4626927103dc697b1084 100644 --- a/MeshLib/Elements/Cell.h +++ b/MeshLib/Elements/Cell.h @@ -38,9 +38,6 @@ public: /// Get the volume of this 3d element. virtual double getVolume() const { return _volume; }; - /// Returns true if the cell is somewhere on the mesh surface and false otherwise. - bool isOnSurface() const; - /// Destructor virtual ~Cell(); diff --git a/MeshLib/MeshSurfaceExtraction.cpp b/MeshLib/MeshSurfaceExtraction.cpp index 26590904b2ec81f9e3c2d5a1ba3339725ca5c8c6..16e1944642fcd516ef2be328b9bb00722ef28c5d 100644 --- a/MeshLib/MeshSurfaceExtraction.cpp +++ b/MeshLib/MeshSurfaceExtraction.cpp @@ -1,5 +1,5 @@ /** - * \file + * \file MeshSurfaceExtraction.cpp * \author Karsten Rink * \date 2013-04-04 * \brief Implementation of the MeshSurfaceExtraction class. @@ -14,7 +14,7 @@ #include "MeshSurfaceExtraction.h" -#include <cassert> +#include <boost/math/constants/constants.hpp> #include "logog/include/logog.hpp" @@ -28,73 +28,77 @@ namespace MeshLib { -void MeshSurfaceExtraction::getSurfaceAreaForNodes(const MeshLib::Mesh &mesh, std::vector<double> &node_area_vec) +std::vector<double> MeshSurfaceExtraction::getSurfaceAreaForNodes(const MeshLib::Mesh &mesh) { - if (mesh.getDimension() == 2) + std::vector<double> node_area_vec; + if (mesh.getDimension() != 2) { - double total_area (0); + ERR ("Error in MeshSurfaceExtraction::getSurfaceAreaForNodes() - Given mesh is no surface mesh (dimension != 2)."); + return node_area_vec; + } - // for each node, a vector containing all the element idget every element - const std::vector<MeshLib::Node*> &nodes = mesh.getNodes(); - const size_t nNodes ( mesh.getNNodes() ); - node_area_vec.reserve(nNodes); - for (size_t n=0; n<nNodes; ++n) - { - double node_area (0); + double total_area (0); - std::vector<MeshLib::Element*> conn_elems = nodes[n]->getElements(); - const size_t nConnElems (conn_elems.size()); + // for each node, a vector containing all the element idget every element + const std::vector<MeshLib::Node*> &nodes = mesh.getNodes(); + const std::size_t nNodes ( mesh.getNNodes() ); + for (std::size_t n=0; n<nNodes; ++n) + { + double node_area (0); - for (size_t i=0; i<nConnElems; ++i) - { - const MeshLib::Element* elem (conn_elems[i]); - const unsigned nElemParts = (elem->getGeomType() == MeshElemType::TRIANGLE) ? 3 : 4; - const double area = conn_elems[i]->getContent() / nElemParts; - node_area += area; - total_area += area; - } + std::vector<MeshLib::Element*> conn_elems = nodes[n]->getElements(); + const std::size_t nConnElems (conn_elems.size()); - node_area_vec.push_back(node_area); + for (std::size_t i=0; i<nConnElems; ++i) + { + const MeshLib::Element* elem (conn_elems[i]); + const unsigned nElemParts = (elem->getGeomType() == MeshElemType::TRIANGLE) ? 3 : 4; + const double area = conn_elems[i]->getContent() / nElemParts; + node_area += area; + total_area += area; } - INFO ("Total surface Area: %f", total_area); + node_area_vec.push_back(node_area); } - else - ERR ("Error in MeshSurfaceExtraction::getSurfaceAreaForNodes() - Given mesh is no surface mesh (dimension != 2)."); + + INFO ("Total surface Area: %f", total_area); + + return node_area_vec; } -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) { - INFO ("Extracting mesh surface..."); - const std::vector<MeshLib::Element*> all_elements (mesh.getElements()); - const std::vector<MeshLib::Node*> all_nodes (mesh.getNodes()); + if (angle< 0 || angle > 90) + { + ERR ("Supported angle between 0 and 90 degrees only."); + return nullptr; + } + INFO ("Extracting mesh surface..."); std::vector<MeshLib::Element*> sfc_elements; - get2DSurfaceElements(all_elements, sfc_elements, dir, mesh.getDimension()); + get2DSurfaceElements(mesh.getElements(), sfc_elements, dir, angle, mesh.getDimension()); if (sfc_elements.empty()) return nullptr; std::vector<MeshLib::Node*> sfc_nodes; - std::vector<unsigned> node_id_map(mesh.getNNodes()); - get2DSurfaceNodes(all_nodes, sfc_nodes, sfc_elements, node_id_map); + std::vector<std::size_t> node_id_map(mesh.getNNodes()); + get2DSurfaceNodes(sfc_nodes, mesh.getNNodes(), sfc_elements, node_id_map); // create new elements vector with newly created nodes std::vector<MeshLib::Element*> new_elements; new_elements.reserve(sfc_elements.size()); - for (auto elem = sfc_elements.begin(); elem != sfc_elements.end(); ++elem) + for (auto elem = sfc_elements.cbegin(); elem != sfc_elements.cend(); ++elem) { - if ((*elem)->getGeomType() == MeshElemType::TRIANGLE) { - MeshLib::Node** tri_nodes = new MeshLib::Node*[3]; - for (unsigned k(0); k<3; k++) - tri_nodes[k] = sfc_nodes[node_id_map[(*elem)->getNode(k)->getID()]]; - new_elements.push_back(new MeshLib::Tri(tri_nodes)); - } else { + unsigned const n_elem_nodes ((*elem)->getNNodes()); + MeshLib::Node** new_nodes = new MeshLib::Node*[n_elem_nodes]; + for (unsigned k(0); k<n_elem_nodes; k++) + new_nodes[k] = sfc_nodes[node_id_map[(*elem)->getNode(k)->getID()]]; + if ((*elem)->getGeomType() == MeshElemType::TRIANGLE) + new_elements.push_back(new MeshLib::Tri(new_nodes)); + else { assert((*elem)->getGeomType() == MeshElemType::QUAD); - MeshLib::Node** quad_nodes = new MeshLib::Node*[4]; - for (unsigned k(0); k<4; k++) - quad_nodes[k] = sfc_nodes[node_id_map[(*elem)->getNode(k)->getID()]]; - new_elements.push_back(new MeshLib::Quad(quad_nodes)); + new_elements.push_back(new MeshLib::Quad(new_nodes)); } delete *elem; } @@ -108,22 +112,24 @@ MeshLib::Mesh* MeshSurfaceExtraction::getMeshSurface(const MeshLib::Mesh &mesh, } MeshLib::Mesh* result (new Mesh(mesh.getName()+"-Surface", sfc_nodes, new_elements)); if (keepOriginalNodeIds) - for (auto node = sfc_nodes.begin(); node != sfc_nodes.end(); ++node) + for (auto node = sfc_nodes.cbegin(); node != sfc_nodes.cend(); ++node) (*node)->setID(id_map[(*node)->getID()]); 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); - bool complete_surface (true); - if (MathLib::scalarProduct(dir, dir) != 0) - complete_surface = false; + bool const complete_surface = (MathLib::scalarProduct(dir, dir) == 0); - for (auto elem = all_elements.begin(); elem != all_elements.end(); ++elem) + double const pi (boost::math::constants::pi<double>()); + double const cos_theta (std::cos(angle * pi / 180.0)); + MathLib::Vector3 const norm_dir (dir.getNormalizedVector()); + + for (auto elem = all_elements.cbegin(); elem != all_elements.cend(); ++elem) { const unsigned element_dimension ((*elem)->getDimension()); if (element_dimension < mesh_dimension) @@ -133,31 +139,31 @@ 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) + MeshLib::Face* face = static_cast<MeshLib::Face*>(*elem); + if (MathLib::scalarProduct(face->getSurfaceNormal().getNormalizedVector(), norm_dir) > cos_theta) continue; } sfc_elements.push_back(*elem); } else { - const MeshLib::Cell* cell = static_cast<MeshLib::Cell*>(*elem); - if (!cell->isOnSurface()) + if (!(*elem)->isBoundaryElement()) continue; - const unsigned nFaces (cell->getNFaces()); + const unsigned nFaces ((*elem)->getNFaces()); for (unsigned j=0; j<nFaces; ++j) { - if (cell->getNeighbor(j) != nullptr) + if ((*elem)->getNeighbor(j) != nullptr) continue; - const MeshLib::Face* face = static_cast<const MeshLib::Face*>(cell->getFace(j)); + const MeshLib::Face* face = static_cast<const MeshLib::Face*>((*elem)->getFace(j)); if (!complete_surface) - if (MathLib::scalarProduct(face->getSurfaceNormal(), dir) <= 0) + { + if (MathLib::scalarProduct(face->getSurfaceNormal().getNormalizedVector(), 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 @@ -167,12 +173,11 @@ void MeshSurfaceExtraction::get2DSurfaceElements(const std::vector<MeshLib::Elem } } -void MeshSurfaceExtraction::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) +void MeshSurfaceExtraction::get2DSurfaceNodes(std::vector<MeshLib::Node*> &sfc_nodes, std::size_t n_all_nodes, const std::vector<MeshLib::Element*> &sfc_elements, std::vector<std::size_t> &node_id_map) { - const size_t nNewElements (sfc_elements.size()); - std::vector<const MeshLib::Node*> tmp_nodes(all_nodes.size(), NULL); - const size_t nNodes (tmp_nodes.size()); - for (unsigned i=0; i<nNewElements; ++i) + const std::size_t nNewElements (sfc_elements.size()); + std::vector<const MeshLib::Node*> tmp_nodes(n_all_nodes, nullptr); + for (std::size_t i=0; i<nNewElements; ++i) { const MeshLib::Element* elem (sfc_elements[i]); for (unsigned j=0; j<elem->getNNodes(); ++j) @@ -181,6 +186,7 @@ void MeshSurfaceExtraction::get2DSurfaceNodes(const std::vector<MeshLib::Node*> tmp_nodes[node->getID()] = node; } } + const std::size_t nNodes (tmp_nodes.size()); for (unsigned i=0; i<nNodes; ++i) { if (tmp_nodes[i]) @@ -191,26 +197,22 @@ 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(mesh.getElements(), sfc_elements, dir, angle, mesh.getDimension()); std::vector<MeshLib::Node*> sfc_nodes; - std::vector<unsigned> node_id_map(mesh.getNNodes()); - get2DSurfaceNodes(all_nodes, sfc_nodes, sfc_elements, node_id_map); + std::vector<std::size_t> node_id_map(mesh.getNNodes()); + get2DSurfaceNodes(sfc_nodes, mesh.getNNodes(), sfc_elements, node_id_map); - const unsigned nElements (sfc_elements.size()); - for (unsigned i=0; i<nElements; ++i) - delete sfc_elements[i]; + for (auto e : sfc_elements) + delete e; - const size_t nNodes (sfc_nodes.size()); + const std::size_t nNodes (sfc_nodes.size()); std::vector<GeoLib::PointWithID*> surface_pnts(nNodes); - for (unsigned i=0; i<nNodes; ++i) + for (std::size_t i=0; i<nNodes; ++i) { surface_pnts[i] = new GeoLib::PointWithID(sfc_nodes[i]->getCoords(), sfc_nodes[i]->getID()); delete sfc_nodes[i]; diff --git a/MeshLib/MeshSurfaceExtraction.h b/MeshLib/MeshSurfaceExtraction.h index ff528d12d09b052fc9cd32d23b1f16e507f6d519..fd58ff99a6738951d658f7212afb5ea9a24eecd6 100644 --- a/MeshLib/MeshSurfaceExtraction.h +++ b/MeshLib/MeshSurfaceExtraction.h @@ -1,5 +1,5 @@ /** - * \file + * \file MeshSurfaceExtraction.h * \author Karsten Rink * \date 2013-04-04 * \brief Definition of the MeshSurfaceExtraction class @@ -36,27 +36,28 @@ class Node; class MeshSurfaceExtraction { public: - /// Returns the area assigned to each node on a surface mesh. - static void getSurfaceAreaForNodes(const MeshLib::Mesh &mesh, std::vector<double> &node_area_vec); + /// Returns a vector of the areas assigned to each node on a surface mesh. + static std::vector<double> getSurfaceAreaForNodes(const MeshLib::Mesh &mesh); /// 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); + static void get2DSurfaceNodes(std::vector<MeshLib::Node*> &sfc_nodes, std::size_t n_all_nodes, const std::vector<MeshLib::Element*> &sfc_elements, std::vector<std::size_t> &node_id_map); }; } // end namespace MeshLib diff --git a/MeshLib/Node.h b/MeshLib/Node.h index 557d26bbdccc7614eebfc80fd28b0e4243c5829f..927bc0c17e11dc602bb17acf09746b61787a210a 100644 --- a/MeshLib/Node.h +++ b/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;