diff --git a/MeshLib/CoordinateSystem.cpp b/MeshLib/CoordinateSystem.cpp index 2cf10443648eed0a62c475d8cef54adec5da674a..f50feec578f1b11459e04253eface2949041c4a7 100644 --- a/MeshLib/CoordinateSystem.cpp +++ b/MeshLib/CoordinateSystem.cpp @@ -9,19 +9,22 @@ #include "CoordinateSystem.h" -#include "MeshLib/Node.h" #include "MeshLib/Elements/Element.h" +#include "MeshLib/Node.h" namespace MeshLib { - -CoordinateSystem::CoordinateSystem(const Element &ele) +CoordinateSystem::CoordinateSystem(const Element& ele) { - GeoLib::AABB const aabb(ele.getNodes(), ele.getNodes() + ele.getNumberOfNodes()); + GeoLib::AABB const aabb(ele.getNodes(), + ele.getNodes() + ele.getNumberOfNodes()); CoordinateSystem const bboxCoordSys(getCoordinateSystem(aabb)); - if (bboxCoordSys.getDimension() >= ele.getDimension()) { + if (bboxCoordSys.getDimension() >= ele.getDimension()) + { _type = bboxCoordSys.getType(); - } else { // e.g. zero volume elements + } + else + { // e.g. zero volume elements if (ele.getDimension() >= 1) { _type = CoordinateSystemType::X; @@ -37,7 +40,8 @@ CoordinateSystem::CoordinateSystem(const Element &ele) } } -unsigned char CoordinateSystem::getCoordinateSystem(const GeoLib::AABB &bbox) const +unsigned char CoordinateSystem::getCoordinateSystem( + const GeoLib::AABB& bbox) const { unsigned char coords = 0; diff --git a/MeshLib/ElementCoordinatesMappingLocal.cpp b/MeshLib/ElementCoordinatesMappingLocal.cpp index 730cdfe732dddf4c3763a0ef857e3fc5e4fc4459..34692d62cae34d5191a76ddd56b6b519296b536e 100644 --- a/MeshLib/ElementCoordinatesMappingLocal.cpp +++ b/MeshLib/ElementCoordinatesMappingLocal.cpp @@ -9,15 +9,14 @@ #include "ElementCoordinatesMappingLocal.h" -#include <limits> #include <cassert> +#include <limits> #include "GeoLib/AnalyticalGeometry.h" - -#include "MeshLib/Elements/Element.h" -#include "MeshLib/Node.h" #include "MathLib/MathTools.h" #include "MathLib/Point3d.h" +#include "MeshLib/Elements/Element.h" +#include "MeshLib/Node.h" namespace detail { diff --git a/MeshLib/ElementStatus.cpp b/MeshLib/ElementStatus.cpp index 5747a9851196f649ca76669b66768546e7a063b8..b82094acd83f99a1c5c3b27a1615075cd46904ed 100644 --- a/MeshLib/ElementStatus.cpp +++ b/MeshLib/ElementStatus.cpp @@ -14,12 +14,12 @@ #include "ElementStatus.h" +#include "Elements/Element.h" #include "Mesh.h" #include "MeshLib/Node.h" -#include "Elements/Element.h" - -namespace MeshLib { +namespace MeshLib +{ ElementStatus::ElementStatus(Mesh const* const mesh, bool hasAnyInactive) : _mesh(mesh), _element_status(mesh->getNumberOfElements(), true), @@ -39,9 +39,12 @@ ElementStatus::ElementStatus(Mesh const* const mesh, { auto* const materialIds = mesh->getProperties().getPropertyVector<int>("MaterialIDs"); - for (auto material_id : vec_inactive_matIDs) { - for (auto e : _mesh->getElements()) { - if ((*materialIds)[e->getID()] == material_id) { + for (auto material_id : vec_inactive_matIDs) + { + for (auto e : _mesh->getElements()) + { + if ((*materialIds)[e->getID()] == material_id) + { setElementStatus(e->getID(), false); } } @@ -49,7 +52,7 @@ ElementStatus::ElementStatus(Mesh const* const mesh, } _vec_active_eles.reserve(getNumberOfActiveElements()); - const std::size_t nElems (_mesh->getNumberOfElements()); + const std::size_t nElems(_mesh->getNumberOfElements()); for (std::size_t i = 0; i < nElems; ++i) { if (_element_status[i]) @@ -60,7 +63,7 @@ ElementStatus::ElementStatus(Mesh const* const mesh, } _vec_active_nodes.reserve(this->getNumberOfActiveNodes()); - const std::size_t nNodes (_mesh->getNumberOfNodes()); + const std::size_t nNodes(_mesh->getNumberOfNodes()); for (std::size_t i = 0; i < nNodes; ++i) { if (_active_nodes[i] > 0) @@ -100,7 +103,8 @@ std::vector<MeshLib::Node*> const& ElementStatus::getActiveNodes() const std::size_t ElementStatus::getNumberOfActiveNodes() const { - return _active_nodes.size() - std::count(_active_nodes.cbegin(), _active_nodes.cend(), 0); + return _active_nodes.size() - + std::count(_active_nodes.cbegin(), _active_nodes.cend(), 0); } std::size_t ElementStatus::getNumberOfActiveElements() const @@ -115,9 +119,10 @@ void ElementStatus::setElementStatus(std::size_t i, bool status) { const int change = (status) ? 1 : -1; _element_status[i] = status; - const unsigned nElemNodes (_mesh->getElement(i)->getNumberOfNodes()); - MeshLib::Node const*const*const nodes = _mesh->getElement(i)->getNodes(); - for (unsigned j=0; j<nElemNodes; ++j) + const unsigned nElemNodes(_mesh->getElement(i)->getNumberOfNodes()); + MeshLib::Node const* const* const nodes = + _mesh->getElement(i)->getNodes(); + for (unsigned j = 0; j < nElemNodes; ++j) { assert(_active_nodes[j] < 255); // if one node has >255 connected // elements the data type is too @@ -129,7 +134,7 @@ void ElementStatus::setElementStatus(std::size_t i, bool status) bool ElementStatus::isActiveNode(MeshLib::Node const* node) const { - return _active_nodes[node->getID()]>0; + return _active_nodes[node->getID()] > 0; } } // namespace MeshLib diff --git a/MeshLib/Elements/CellRule.cpp b/MeshLib/Elements/CellRule.cpp index bce43895fee919c95ef8ab7e08af7518b3279783..17e0a3fd4b17d88180b1f59bf03ba08a383625e7 100644 --- a/MeshLib/Elements/CellRule.cpp +++ b/MeshLib/Elements/CellRule.cpp @@ -10,12 +10,12 @@ #include "CellRule.h" -#include "MeshLib/Node.h" #include "Element.h" #include "FaceRule.h" +#include "MeshLib/Node.h" -namespace MeshLib { - +namespace MeshLib +{ bool CellRule::testElementNodeOrder(const Element* e) { Eigen::Vector3d const cc = @@ -24,10 +24,10 @@ bool CellRule::testElementNodeOrder(const Element* e) for (unsigned j = 0; j < nFaces; ++j) { MeshLib::Element const* const face(e->getFace(j)); - // Node 1 is checked below because that way all nodes are used for the test - // at some point, while for node 0 at least one node in every element - // type would be used for checking twice and one wouldn't be checked at - // all. (based on the definition of the _face_nodes variable) + // Node 1 is checked below because that way all nodes are used for the + // test at some point, while for node 0 at least one node in every + // element type would be used for checking twice and one wouldn't be + // checked at all. (based on the definition of the _face_nodes variable) auto const x = Eigen::Map<Eigen::Vector3d const>(face->getNode(1)->getCoords()); Eigen::Vector3d const cx = x - cc; diff --git a/MeshLib/Elements/EdgeReturn.cpp b/MeshLib/Elements/EdgeReturn.cpp index fa94ff5d60910c9628ac20c2e9bad5c1dbf97b81..500ef0e14e25f4322facc11dddf8a548c2020183 100644 --- a/MeshLib/Elements/EdgeReturn.cpp +++ b/MeshLib/Elements/EdgeReturn.cpp @@ -11,21 +11,19 @@ #include "EdgeReturn.h" #include "BaseLib/Logging.h" - -#include "MeshLib/Node.h" #include "Element.h" #include "Line.h" +#include "MeshLib/Node.h" namespace MeshLib { - const Element* LinearEdgeReturn::getEdge(const Element* e, unsigned i) { if (i < e->getNumberOfEdges()) { auto** nodes = new Node*[2]; - nodes[0] = const_cast<Node*>(e->getEdgeNode(i,0)); - nodes[1] = const_cast<Node*>(e->getEdgeNode(i,1)); + nodes[0] = const_cast<Node*>(e->getEdgeNode(i, 0)); + nodes[1] = const_cast<Node*>(e->getEdgeNode(i, 1)); return new Line(nodes, e->getID()); } ERR("Error in MeshLib::Element::getEdge() - Index does not exist."); @@ -37,9 +35,9 @@ const Element* QuadraticEdgeReturn::getEdge(const Element* e, unsigned i) if (i < e->getNumberOfEdges()) { auto** nodes = new Node*[3]; - nodes[0] = const_cast<Node*>(e->getEdgeNode(i,0)); - nodes[1] = const_cast<Node*>(e->getEdgeNode(i,1)); - nodes[2] = const_cast<Node*>(e->getEdgeNode(i,2)); + nodes[0] = const_cast<Node*>(e->getEdgeNode(i, 0)); + nodes[1] = const_cast<Node*>(e->getEdgeNode(i, 1)); + nodes[2] = const_cast<Node*>(e->getEdgeNode(i, 2)); return new Line3(nodes, e->getID()); } ERR("Error in MeshLib::Element::getEdge() - Index does not exist."); diff --git a/MeshLib/Elements/Element.cpp b/MeshLib/Elements/Element.cpp index f73cfa6ce9523a3af2d071f1b84c3e211177441c..95490314afb8d29773f4cf62b010608c46f7b2e6 100644 --- a/MeshLib/Elements/Element.cpp +++ b/MeshLib/Elements/Element.cpp @@ -15,15 +15,13 @@ #include "Element.h" #include "BaseLib/Logging.h" - +#include "Line.h" #include "MathLib/GeometricBasics.h" #include "MathLib/MathTools.h" #include "MeshLib/Node.h" -#include "Line.h" - -namespace MeshLib { - +namespace MeshLib +{ Element::Element(std::size_t id) : _nodes(nullptr), _id(id), _content(-1.0), _neighbors(nullptr) { @@ -31,8 +29,8 @@ Element::Element(std::size_t id) Element::~Element() { - delete [] this->_nodes; - delete [] this->_neighbors; + delete[] this->_nodes; + delete[] this->_neighbors; } void Element::setNeighbor(Element* neighbor, unsigned const face_id) @@ -58,11 +56,11 @@ std::optional<unsigned> Element::addNeighbor(Element* e) } Node* face_nodes[3]; - const unsigned nNodes (this->getNumberOfBaseNodes()); - const unsigned eNodes (e->getNumberOfBaseNodes()); + const unsigned nNodes(this->getNumberOfBaseNodes()); + const unsigned eNodes(e->getNumberOfBaseNodes()); const Node* const* e_nodes = e->getNodes(); unsigned count(0); - const unsigned dim (this->getDimension()); + const unsigned dim(this->getDimension()); for (unsigned i(0); i < nNodes; i++) { for (unsigned j(0); j < eNodes; j++) @@ -70,10 +68,11 @@ std::optional<unsigned> Element::addNeighbor(Element* e) if (_nodes[i] == e_nodes[j]) { face_nodes[count] = _nodes[i]; - // increment shared nodes counter and check if enough nodes are similar to be sure e is a neighbour of this - if ((++count)>=dim) + // increment shared nodes counter and check if enough nodes are + // similar to be sure e is a neighbour of this + if ((++count) >= dim) { - _neighbors[ this->identifyFace(face_nodes) ] = e; + _neighbors[this->identifyFace(face_nodes)] = e; return std::optional<unsigned>(e->identifyFace(face_nodes)); } } @@ -99,7 +98,7 @@ const Element* Element::getNeighbor(unsigned i) const unsigned Element::getNodeIDinElement(const MeshLib::Node* node) const { - const unsigned nNodes (this->getNumberOfNodes()); + const unsigned nNodes(this->getNumberOfNodes()); for (unsigned i(0); i < nNodes; i++) { if (node == _nodes[i]) @@ -151,21 +150,22 @@ std::size_t Element::getNodeIndex(unsigned i) const bool Element::isBoundaryElement() const { - return std::any_of(_neighbors, _neighbors + this->getNumberOfNeighbors(), - [](MeshLib::Element const*const e){ return e == nullptr; }); + return std::any_of( + _neighbors, _neighbors + this->getNumberOfNeighbors(), + [](MeshLib::Element const* const e) { return e == nullptr; }); } #ifndef NDEBUG std::ostream& operator<<(std::ostream& os, Element const& e) { - os << "Element #" << e._id << " @ " << &e << " with " << e.getNumberOfNeighbors() - << " neighbours\n"; + os << "Element #" << e._id << " @ " << &e << " with " + << e.getNumberOfNeighbors() << " neighbours\n"; unsigned const nnodes = e.getNumberOfNodes(); MeshLib::Node* const* const nodes = e.getNodes(); os << "MeshElemType: " - << static_cast<std::underlying_type<MeshElemType>::type>(e.getGeomType()) - << " with " << nnodes << " nodes: { "; + << static_cast<std::underlying_type<MeshElemType>::type>(e.getGeomType()) + << " with " << nnodes << " nodes: { "; for (unsigned n = 0; n < nnodes; ++n) { os << nodes[n]->getID() << " @ " << nodes[n] << " "; @@ -246,7 +246,8 @@ std::pair<double, double> computeSqrEdgeLengthRange(Element const& element) bool isPointInElementXY(MathLib::Point3d const& p, Element const& e) { - for(std::size_t i(0); i<e.getNumberOfBaseNodes(); ++i) { + for (std::size_t i(0); i < e.getNumberOfBaseNodes(); ++i) + { if (MathLib::sqrDist2d(p, *e.getNode(i)) < std::numeric_limits<double>::epsilon()) { diff --git a/MeshLib/Elements/FaceRule.cpp b/MeshLib/Elements/FaceRule.cpp index 5a56c3902df69c26751177e38a21eb021ff69542..7471c4a3457d72bef27414281a3447601ca4290b 100644 --- a/MeshLib/Elements/FaceRule.cpp +++ b/MeshLib/Elements/FaceRule.cpp @@ -10,9 +10,9 @@ #include "FaceRule.h" +#include "Element.h" #include "MathLib/MathTools.h" #include "MeshLib/Node.h" -#include "Element.h" namespace MeshLib { diff --git a/MeshLib/Elements/HexRule20.cpp b/MeshLib/Elements/HexRule20.cpp index 71bbd5625921374e1c4c3a00d12a554581edf260..0f815d33779e6acb2aca1c81296f433e6ba144f8 100644 --- a/MeshLib/Elements/HexRule20.cpp +++ b/MeshLib/Elements/HexRule20.cpp @@ -13,12 +13,12 @@ #include <array> #include "BaseLib/Logging.h" - +#include "Line.h" #include "MeshLib/Node.h" #include "Quad.h" -#include "Line.h" -namespace MeshLib { +namespace MeshLib +{ const unsigned HexRule20::face_nodes[6][8] = { {0, 3, 2, 1, 11, 10, 9, 8}, // Face 0 {0, 1, 5, 4, 8, 17, 12, 16}, // Face 1 @@ -28,20 +28,19 @@ const unsigned HexRule20::face_nodes[6][8] = { {4, 5, 6, 7, 12, 13, 14, 15} // Face 5 }; -const unsigned HexRule20::edge_nodes[12][3] = -{ - {0, 1, 8}, // Edge 0 - {1, 2, 9}, // Edge 1 - {2, 3, 10}, // Edge 2 - {0, 3, 11}, // Edge 3 - {4, 5, 12}, // Edge 4 - {5, 6, 13}, // Edge 5 - {6, 7, 14}, // Edge 6 - {4, 7, 15}, // Edge 7 - {0, 4, 16}, // Edge 8 - {1, 5, 17}, // Edge 9 - {2, 6, 18}, // Edge 10 - {3, 7, 19} // Edge 11 +const unsigned HexRule20::edge_nodes[12][3] = { + {0, 1, 8}, // Edge 0 + {1, 2, 9}, // Edge 1 + {2, 3, 10}, // Edge 2 + {0, 3, 11}, // Edge 3 + {4, 5, 12}, // Edge 4 + {5, 6, 13}, // Edge 5 + {6, 7, 14}, // Edge 6 + {4, 7, 15}, // Edge 7 + {0, 4, 16}, // Edge 8 + {1, 5, 17}, // Edge 9 + {2, 6, 18}, // Edge 10 + {3, 7, 19} // Edge 11 }; const Element* HexRule20::getFace(const Element* e, unsigned i) @@ -59,4 +58,4 @@ const Element* HexRule20::getFace(const Element* e, unsigned i) return nullptr; } -} // end namespace MeshLib +} // end namespace MeshLib diff --git a/MeshLib/Elements/HexRule8.cpp b/MeshLib/Elements/HexRule8.cpp index 63bb68b7dbaca8ae56288ad515f3a8cdf4f0df95..394bcabeca1a929bf3ba51d13fcc823ed31d203e 100644 --- a/MeshLib/Elements/HexRule8.cpp +++ b/MeshLib/Elements/HexRule8.cpp @@ -13,39 +13,35 @@ #include <array> #include "BaseLib/Logging.h" - +#include "Line.h" #include "MathLib/GeometricBasics.h" - #include "MeshLib/Node.h" #include "Quad.h" -#include "Line.h" - -namespace MeshLib { -const unsigned HexRule8::face_nodes[6][4] = +namespace MeshLib { - {0, 3, 2, 1}, // Face 0 - {0, 1, 5, 4}, // Face 1 - {1, 2, 6, 5}, // Face 2 - {2, 3, 7, 6}, // Face 3 - {3, 0, 4, 7}, // Face 4 - {4, 5, 6, 7} // Face 5 +const unsigned HexRule8::face_nodes[6][4] = { + {0, 3, 2, 1}, // Face 0 + {0, 1, 5, 4}, // Face 1 + {1, 2, 6, 5}, // Face 2 + {2, 3, 7, 6}, // Face 3 + {3, 0, 4, 7}, // Face 4 + {4, 5, 6, 7} // Face 5 }; -const unsigned HexRule8::edge_nodes[12][2] = -{ - {0, 1}, // Edge 0 - {1, 2}, // Edge 1 - {2, 3}, // Edge 2 - {0, 3}, // Edge 3 - {4, 5}, // Edge 4 - {5, 6}, // Edge 5 - {6, 7}, // Edge 6 - {4, 7}, // Edge 7 - {0, 4}, // Edge 8 - {1, 5}, // Edge 9 - {2, 6}, // Edge 10 - {3, 7} // Edge 11 +const unsigned HexRule8::edge_nodes[12][2] = { + {0, 1}, // Edge 0 + {1, 2}, // Edge 1 + {2, 3}, // Edge 2 + {0, 3}, // Edge 3 + {4, 5}, // Edge 4 + {5, 6}, // Edge 5 + {6, 7}, // Edge 6 + {4, 7}, // Edge 7 + {0, 4}, // Edge 8 + {1, 5}, // Edge 9 + {2, 6}, // Edge 10 + {3, 7} // Edge 11 }; const Element* HexRule8::getFace(const Element* e, unsigned i) @@ -65,12 +61,18 @@ const Element* HexRule8::getFace(const Element* e, unsigned i) double HexRule8::computeVolume(Node const* const* _nodes) { - return MathLib::calcTetrahedronVolume(*_nodes[4], *_nodes[7], *_nodes[5], *_nodes[0]) - + MathLib::calcTetrahedronVolume(*_nodes[5], *_nodes[3], *_nodes[1], *_nodes[0]) - + MathLib::calcTetrahedronVolume(*_nodes[5], *_nodes[7], *_nodes[3], *_nodes[0]) - + MathLib::calcTetrahedronVolume(*_nodes[5], *_nodes[7], *_nodes[6], *_nodes[2]) - + MathLib::calcTetrahedronVolume(*_nodes[1], *_nodes[3], *_nodes[5], *_nodes[2]) - + MathLib::calcTetrahedronVolume(*_nodes[3], *_nodes[7], *_nodes[5], *_nodes[2]); + return MathLib::calcTetrahedronVolume( + *_nodes[4], *_nodes[7], *_nodes[5], *_nodes[0]) + + MathLib::calcTetrahedronVolume( + *_nodes[5], *_nodes[3], *_nodes[1], *_nodes[0]) + + MathLib::calcTetrahedronVolume( + *_nodes[5], *_nodes[7], *_nodes[3], *_nodes[0]) + + MathLib::calcTetrahedronVolume( + *_nodes[5], *_nodes[7], *_nodes[6], *_nodes[2]) + + MathLib::calcTetrahedronVolume( + *_nodes[1], *_nodes[3], *_nodes[5], *_nodes[2]) + + MathLib::calcTetrahedronVolume( + *_nodes[3], *_nodes[7], *_nodes[5], *_nodes[2]); } bool HexRule8::isPntInElement(Node const* const* nodes, @@ -93,7 +95,7 @@ bool HexRule8::isPntInElement(Node const* const* nodes, unsigned HexRule8::identifyFace(Node const* const* _nodes, Node* nodes[3]) { - for (unsigned i=0; i<6; i++) + for (unsigned i = 0; i < 6; i++) { unsigned flag(0); for (unsigned j = 0; j < 4; j++) @@ -119,19 +121,19 @@ ElementErrorCode HexRule8::validate(const Element* e) ElementErrorCode error_code; error_code[ElementErrorFlag::ZeroVolume] = hasZeroVolume(*e); - for (unsigned i=0; i<6; ++i) + for (unsigned i = 0; i < 6; ++i) { if (error_code.all()) { break; } - const MeshLib::Element* quad (e->getFace(i)); + const MeshLib::Element* quad(e->getFace(i)); error_code |= quad->validate(); delete quad; } - error_code[ElementErrorFlag::NodeOrder] = !e->testElementNodeOrder(); + error_code[ElementErrorFlag::NodeOrder] = !e->testElementNodeOrder(); return error_code; } -} // end namespace MeshLib +} // end namespace MeshLib diff --git a/MeshLib/Elements/LineRule2.cpp b/MeshLib/Elements/LineRule2.cpp index 6f3fa78e2e09dd0916c450237907ffff37e5ad55..5fc4ad6588a762e16bbb9b0a2be6de80e6330036 100644 --- a/MeshLib/Elements/LineRule2.cpp +++ b/MeshLib/Elements/LineRule2.cpp @@ -13,11 +13,10 @@ #include "MathLib/MathTools.h" #include "MeshLib/Node.h" -namespace MeshLib { - -const unsigned LineRule2::edge_nodes[1][2] = +namespace MeshLib { - {0, 1} // Edge 0 +const unsigned LineRule2::edge_nodes[1][2] = { + {0, 1} // Edge 0 }; double LineRule2::computeVolume(Node const* const* _nodes) @@ -55,4 +54,4 @@ ElementErrorCode LineRule2::validate(const Element* e) return error_code; } -} // end namespace MeshLib +} // end namespace MeshLib diff --git a/MeshLib/Elements/MapBulkElementPoint.cpp b/MeshLib/Elements/MapBulkElementPoint.cpp index a3cfcc73d9c0b8d8ed786dee6cc05a183f5ca72d..d2039ee084fb86c1355c4d72aca10b8d10d6fc09 100644 --- a/MeshLib/Elements/MapBulkElementPoint.cpp +++ b/MeshLib/Elements/MapBulkElementPoint.cpp @@ -8,10 +8,10 @@ * */ -#include <array> - #include "MapBulkElementPoint.h" +#include <array> + namespace MeshLib { MathLib::Point3d getBulkElementPoint(MeshLib::Tri const& /*tri*/, @@ -37,17 +37,20 @@ MathLib::Point3d getBulkElementPoint(MeshLib::Quad const& /*quad*/, std::size_t const face_id, MathLib::WeightedPoint1D const& wp) { - switch (face_id) { - case 0: - return MathLib::Point3d{std::array<double, 3>{{wp[0], 0.0, 0.0}}}; - case 1: - return MathLib::Point3d{std::array<double, 3>{{1.0, wp[0], 0.0}}}; - case 2: - return MathLib::Point3d{std::array<double, 3>{{1.0 - wp[0], 1.0, 0.0}}}; - case 3: - return MathLib::Point3d{std::array<double, 3>{{0.0, 1.0 - wp[0], 0.0}}}; - default: - OGS_FATAL("Invalid face id '{:d}' for the quad.", face_id); + switch (face_id) + { + case 0: + return MathLib::Point3d{std::array<double, 3>{{wp[0], 0.0, 0.0}}}; + case 1: + return MathLib::Point3d{std::array<double, 3>{{1.0, wp[0], 0.0}}}; + case 2: + return MathLib::Point3d{ + std::array<double, 3>{{1.0 - wp[0], 1.0, 0.0}}}; + case 3: + return MathLib::Point3d{ + std::array<double, 3>{{0.0, 1.0 - wp[0], 0.0}}}; + default: + OGS_FATAL("Invalid face id '{:d}' for the quad.", face_id); } } @@ -89,16 +92,15 @@ MathLib::Point3d getBulkElementPoint(MeshLib::Prism const& /*prism*/, std::array<double, 3>{{wp[1], wp[0], -1.0}}}; case 1: return MathLib::Point3d{ - std::array<double, 3>{{wp[0]/2.0 + 0.5, 0.0, wp[1]}}}; + std::array<double, 3>{{wp[0] / 2.0 + 0.5, 0.0, wp[1]}}}; case 2: - return MathLib::Point3d{ - std::array<double, 3>{{0.5 - wp[0]/2.0, 0.5 + wp[0]/2.0, wp[1]}}}; + return MathLib::Point3d{std::array<double, 3>{ + {0.5 - wp[0] / 2.0, 0.5 + wp[0] / 2.0, wp[1]}}}; case 3: return MathLib::Point3d{ - std::array<double, 3>{{0, -wp[0]/2.0 + 0.5, wp[1]}}}; + std::array<double, 3>{{0, -wp[0] / 2.0 + 0.5, wp[1]}}}; case 4: - return MathLib::Point3d{ - std::array<double, 3>{{wp[0], wp[1], 1.0}}}; + return MathLib::Point3d{std::array<double, 3>{{wp[0], wp[1], 1.0}}}; default: OGS_FATAL("Invalid face id '{:d}' for the prism.", face_id); } @@ -120,8 +122,8 @@ MathLib::Point3d getBulkElementPoint(MeshLib::Pyramid const& /*pyramid*/, return MathLib::Point3d{ std::array<double, 3>{{1 - 2 * wp[0], 1.0, 2 * wp[1] - 1}}}; case 3: - return MathLib::Point3d{std::array<double, 3>{ - {-1.0, 2 * wp[1] - 1, 2 * wp[0] - 1}}}; + return MathLib::Point3d{ + std::array<double, 3>{{-1.0, 2 * wp[1] - 1, 2 * wp[0] - 1}}}; case 4: return MathLib::Point3d{ std::array<double, 3>{{-wp[0], wp[1], -1.0}}}; @@ -144,8 +146,7 @@ MathLib::Point3d getBulkElementPoint(MeshLib::Tet const& /*tet*/, return MathLib::Point3d{ std::array<double, 3>{{1 - wp[0] - wp[1], wp[0], wp[1]}}}; case 3: - return MathLib::Point3d{ - std::array<double, 3>{{0, wp[1], wp[0]}}}; + return MathLib::Point3d{std::array<double, 3>{{0, wp[1], wp[0]}}}; default: OGS_FATAL("Invalid face id '{:d}' for the tetrahedron.", face_id); } diff --git a/MeshLib/Elements/PointRule1.cpp b/MeshLib/Elements/PointRule1.cpp index dd13cf1adb5b15695a21b227a21dc7862a4928be..e812fb9372578117147ee8841b5ed986299b9e90 100644 --- a/MeshLib/Elements/PointRule1.cpp +++ b/MeshLib/Elements/PointRule1.cpp @@ -13,12 +13,9 @@ #include "MathLib/Point3d.h" #include "MeshLib/Node.h" -namespace MeshLib { - -const unsigned PointRule1::edge_nodes[1][1] = +namespace MeshLib { - {0} -}; +const unsigned PointRule1::edge_nodes[1][1] = {{0}}; double PointRule1::computeVolume(Node const* const* /*_nodes*/) { @@ -48,4 +45,4 @@ ElementErrorCode PointRule1::validate(const Element* e) return error_code; } -} // end namespace MeshLib +} // end namespace MeshLib diff --git a/MeshLib/Elements/PrismRule15.cpp b/MeshLib/Elements/PrismRule15.cpp index 6ecd50b88045f8576b1698dbf5c037b008f04180..7d1d6c37bb1ccfdb582d0bde54e323ab5ddf5479 100644 --- a/MeshLib/Elements/PrismRule15.cpp +++ b/MeshLib/Elements/PrismRule15.cpp @@ -11,36 +11,33 @@ #include "PrismRule15.h" #include "BaseLib/Logging.h" - #include "MeshLib/Node.h" #include "Quad.h" #include "Tri.h" -namespace MeshLib { - -const unsigned PrismRule15::face_nodes[5][8] = +namespace MeshLib { - {0, 2, 1, 8, 7, 6, 99, 99}, // Face 0 - {0, 1, 4, 3, 6, 10, 12, 9}, // Face 1 - {1, 2, 5, 4, 7, 11, 13, 10}, // Face 2 - {2, 0, 3, 5, 8, 9, 14, 11}, // Face 3 +const unsigned PrismRule15::face_nodes[5][8] = { + {0, 2, 1, 8, 7, 6, 99, 99}, // Face 0 + {0, 1, 4, 3, 6, 10, 12, 9}, // Face 1 + {1, 2, 5, 4, 7, 11, 13, 10}, // Face 2 + {2, 0, 3, 5, 8, 9, 14, 11}, // Face 3 {3, 4, 5, 12, 13, 14, 99, 99} // Face 4 }; -const unsigned PrismRule15::edge_nodes[9][3] = -{ - {0, 1, 6}, // Edge 0 - {1, 2, 7}, // Edge 1 - {0, 2, 8}, // Edge 2 - {0, 3, 9}, // Edge 3 - {1, 4, 10}, // Edge 4 - {2, 5, 11}, // Edge 5 - {3, 4, 12}, // Edge 6 - {4, 5, 13}, // Edge 7 - {3, 5, 14} // Edge 8 +const unsigned PrismRule15::edge_nodes[9][3] = { + {0, 1, 6}, // Edge 0 + {1, 2, 7}, // Edge 1 + {0, 2, 8}, // Edge 2 + {0, 3, 9}, // Edge 3 + {1, 4, 10}, // Edge 4 + {2, 5, 11}, // Edge 5 + {3, 4, 12}, // Edge 6 + {4, 5, 13}, // Edge 7 + {3, 5, 14} // Edge 8 }; -const unsigned PrismRule15::n_face_nodes[5] = { 6, 8, 8, 8, 6 }; +const unsigned PrismRule15::n_face_nodes[5] = {6, 8, 8, 8, 6}; const Element* PrismRule15::getFace(const Element* e, unsigned i) { @@ -64,4 +61,4 @@ const Element* PrismRule15::getFace(const Element* e, unsigned i) return nullptr; } -} // end namespace MeshLib +} // end namespace MeshLib diff --git a/MeshLib/Elements/PrismRule6.cpp b/MeshLib/Elements/PrismRule6.cpp index b46fef3ff361983c31ae3b2338accad9b4eb172e..5bc0de116276cb5913e1012b263637efcf5f0747 100644 --- a/MeshLib/Elements/PrismRule6.cpp +++ b/MeshLib/Elements/PrismRule6.cpp @@ -11,38 +11,34 @@ #include "PrismRule6.h" #include "BaseLib/Logging.h" - #include "MathLib/GeometricBasics.h" - #include "MeshLib/Node.h" #include "Quad.h" #include "Tri.h" -namespace MeshLib { - -const unsigned PrismRule6::face_nodes[5][4] = +namespace MeshLib { - {0, 2, 1, 99}, // Face 0 - {0, 1, 4, 3}, // Face 1 - {1, 2, 5, 4}, // Face 2 - {2, 0, 3, 5}, // Face 3 - {3, 4, 5, 99} // Face 4 +const unsigned PrismRule6::face_nodes[5][4] = { + {0, 2, 1, 99}, // Face 0 + {0, 1, 4, 3}, // Face 1 + {1, 2, 5, 4}, // Face 2 + {2, 0, 3, 5}, // Face 3 + {3, 4, 5, 99} // Face 4 }; -const unsigned PrismRule6::edge_nodes[9][2] = -{ - {0, 1}, // Edge 0 - {1, 2}, // Edge 1 - {0, 2}, // Edge 2 - {0, 3}, // Edge 3 - {1, 4}, // Edge 4 - {2, 5}, // Edge 5 - {3, 4}, // Edge 6 - {4, 5}, // Edge 7 - {3, 5} // Edge 8 +const unsigned PrismRule6::edge_nodes[9][2] = { + {0, 1}, // Edge 0 + {1, 2}, // Edge 1 + {0, 2}, // Edge 2 + {0, 3}, // Edge 3 + {1, 4}, // Edge 4 + {2, 5}, // Edge 5 + {3, 4}, // Edge 6 + {4, 5}, // Edge 7 + {3, 5} // Edge 8 }; -const unsigned PrismRule6::n_face_nodes[5] = { 3, 4, 4, 4, 3 }; +const unsigned PrismRule6::n_face_nodes[5] = {3, 4, 4, 4, 3}; const Element* PrismRule6::getFace(const Element* e, unsigned i) { @@ -68,9 +64,12 @@ const Element* PrismRule6::getFace(const Element* e, unsigned i) double PrismRule6::computeVolume(Node const* const* _nodes) { - return MathLib::calcTetrahedronVolume(*_nodes[0], *_nodes[1], *_nodes[2], *_nodes[3]) - + MathLib::calcTetrahedronVolume(*_nodes[1], *_nodes[4], *_nodes[2], *_nodes[3]) - + MathLib::calcTetrahedronVolume(*_nodes[2], *_nodes[4], *_nodes[5], *_nodes[3]); + return MathLib::calcTetrahedronVolume( + *_nodes[0], *_nodes[1], *_nodes[2], *_nodes[3]) + + MathLib::calcTetrahedronVolume( + *_nodes[1], *_nodes[4], *_nodes[2], *_nodes[3]) + + MathLib::calcTetrahedronVolume( + *_nodes[2], *_nodes[4], *_nodes[5], *_nodes[3]); } bool PrismRule6::isPntInElement(Node const* const* nodes, @@ -87,7 +86,7 @@ bool PrismRule6::isPntInElement(Node const* const* nodes, unsigned PrismRule6::identifyFace(Node const* const* _nodes, Node* nodes[3]) { - for (unsigned i=0; i<5; i++) + for (unsigned i = 0; i < 5; i++) { unsigned flag(0); for (unsigned j = 0; j < 4; j++) @@ -114,7 +113,7 @@ ElementErrorCode PrismRule6::validate(const Element* e) ElementErrorCode error_code; error_code[ElementErrorFlag::ZeroVolume] = hasZeroVolume(*e); - for (unsigned i=1; i<4; ++i) + for (unsigned i = 1; i < 4; ++i) { const auto* quad(dynamic_cast<const MeshLib::Quad*>(e->getFace(i))); if (quad) @@ -131,4 +130,4 @@ ElementErrorCode PrismRule6::validate(const Element* e) return error_code; } -} // end namespace MeshLib +} // end namespace MeshLib diff --git a/MeshLib/Elements/PyramidRule13.cpp b/MeshLib/Elements/PyramidRule13.cpp index 02ecd73baa9e3c8f0f0da22d3640dac12978b3fc..9d76fac3e951fe2d2a7c6780c5e77bdb2f0d4d36 100644 --- a/MeshLib/Elements/PyramidRule13.cpp +++ b/MeshLib/Elements/PyramidRule13.cpp @@ -11,39 +11,36 @@ #include "PyramidRule13.h" #include "BaseLib/Logging.h" - #include "MeshLib/Node.h" #include "Quad.h" #include "Tri.h" -namespace MeshLib { - -const unsigned PyramidRule13::face_nodes[5][8] = +namespace MeshLib { - {0, 1, 4, 5, 10, 9, 99, 99}, // Face 0 - {1, 2, 4, 6, 11, 10, 99, 99}, // Face 1 - {2, 3, 4, 7, 12, 11, 99, 99}, // Face 2 - {3, 0, 4, 8, 9, 12, 99, 99}, // Face 3 - {0, 3, 2, 1, 8, 7, 6, 5} // Face 4 +const unsigned PyramidRule13::face_nodes[5][8] = { + {0, 1, 4, 5, 10, 9, 99, 99}, // Face 0 + {1, 2, 4, 6, 11, 10, 99, 99}, // Face 1 + {2, 3, 4, 7, 12, 11, 99, 99}, // Face 2 + {3, 0, 4, 8, 9, 12, 99, 99}, // Face 3 + {0, 3, 2, 1, 8, 7, 6, 5} // Face 4 }; -const unsigned PyramidRule13::edge_nodes[8][3] = -{ - {0, 1, 5}, // Edge 0 - {1, 2, 6}, // Edge 1 - {2, 3, 7}, // Edge 2 - {0, 3, 8}, // Edge 3 - {0, 4, 9}, // Edge 4 - {1, 4, 10}, // Edge 5 - {2, 4, 11}, // Edge 6 - {3, 4, 12} // Edge 7 +const unsigned PyramidRule13::edge_nodes[8][3] = { + {0, 1, 5}, // Edge 0 + {1, 2, 6}, // Edge 1 + {2, 3, 7}, // Edge 2 + {0, 3, 8}, // Edge 3 + {0, 4, 9}, // Edge 4 + {1, 4, 10}, // Edge 5 + {2, 4, 11}, // Edge 6 + {3, 4, 12} // Edge 7 }; -const unsigned PyramidRule13::n_face_nodes[5] = { 6, 6, 6, 6, 8 }; +const unsigned PyramidRule13::n_face_nodes[5] = {6, 6, 6, 6, 8}; const Element* PyramidRule13::getFace(const Element* e, unsigned i) { - if (i<e->getNumberOfFaces()) + if (i < e->getNumberOfFaces()) { unsigned nFaceNodes(PyramidRule13::n_face_nodes[i]); auto** nodes = new Node*[nFaceNodes]; @@ -63,4 +60,4 @@ const Element* PyramidRule13::getFace(const Element* e, unsigned i) return nullptr; } -} // end namespace MeshLib +} // end namespace MeshLib diff --git a/MeshLib/Elements/PyramidRule5.cpp b/MeshLib/Elements/PyramidRule5.cpp index dbc56555705955337d5d14dffe90690b7d1a6dfc..7b6486da71bf606014816c26fa9e2a7639d3cec9 100644 --- a/MeshLib/Elements/PyramidRule5.cpp +++ b/MeshLib/Elements/PyramidRule5.cpp @@ -11,41 +11,37 @@ #include "PyramidRule5.h" #include "BaseLib/Logging.h" - #include "MathLib/GeometricBasics.h" - #include "MeshLib/Node.h" #include "Quad.h" #include "Tri.h" -namespace MeshLib { - -const unsigned PyramidRule5::face_nodes[5][4] = +namespace MeshLib { - {0, 1, 4, 99}, // Face 0 - {1, 2, 4, 99}, // Face 1 - {2, 3, 4, 99}, // Face 2 - {3, 0, 4, 99}, // Face 3 - {0, 3, 2, 1} // Face 4 +const unsigned PyramidRule5::face_nodes[5][4] = { + {0, 1, 4, 99}, // Face 0 + {1, 2, 4, 99}, // Face 1 + {2, 3, 4, 99}, // Face 2 + {3, 0, 4, 99}, // Face 3 + {0, 3, 2, 1} // Face 4 }; -const unsigned PyramidRule5::edge_nodes[8][2] = -{ - {0, 1}, // Edge 0 - {1, 2}, // Edge 1 - {2, 3}, // Edge 2 - {0, 3}, // Edge 3 - {0, 4}, // Edge 4 - {1, 4}, // Edge 5 - {2, 4}, // Edge 6 - {3, 4} // Edge 7 +const unsigned PyramidRule5::edge_nodes[8][2] = { + {0, 1}, // Edge 0 + {1, 2}, // Edge 1 + {2, 3}, // Edge 2 + {0, 3}, // Edge 3 + {0, 4}, // Edge 4 + {1, 4}, // Edge 5 + {2, 4}, // Edge 6 + {3, 4} // Edge 7 }; -const unsigned PyramidRule5::n_face_nodes[5] = { 3, 3, 3, 3, 4 }; +const unsigned PyramidRule5::n_face_nodes[5] = {3, 3, 3, 3, 4}; const Element* PyramidRule5::getFace(const Element* e, unsigned i) { - if (i<e->getNumberOfFaces()) + if (i < e->getNumberOfFaces()) { unsigned nFaceNodes(PyramidRule5::n_face_nodes[i]); auto** nodes = new Node*[nFaceNodes]; @@ -67,8 +63,10 @@ const Element* PyramidRule5::getFace(const Element* e, unsigned i) double PyramidRule5::computeVolume(Node const* const* _nodes) { - return MathLib::calcTetrahedronVolume(*_nodes[0], *_nodes[1], *_nodes[2], *_nodes[4]) - + MathLib::calcTetrahedronVolume(*_nodes[2], *_nodes[3], *_nodes[0], *_nodes[4]); + return MathLib::calcTetrahedronVolume( + *_nodes[0], *_nodes[1], *_nodes[2], *_nodes[4]) + + MathLib::calcTetrahedronVolume( + *_nodes[2], *_nodes[3], *_nodes[0], *_nodes[4]); } bool PyramidRule5::isPntInElement(Node const* const* nodes, @@ -83,7 +81,7 @@ bool PyramidRule5::isPntInElement(Node const* const* nodes, unsigned PyramidRule5::identifyFace(Node const* const* _nodes, Node* nodes[3]) { - for (unsigned i=0; i<5; i++) + for (unsigned i = 0; i < 5; i++) { unsigned flag(0); for (unsigned j = 0; j < 4; j++) @@ -125,4 +123,4 @@ ElementErrorCode PyramidRule5::validate(const Element* e) return error_code; } -} // end namespace MeshLib +} // end namespace MeshLib diff --git a/MeshLib/Elements/QuadRule4.cpp b/MeshLib/Elements/QuadRule4.cpp index b98b106435871af36293b6a558429c061f267ea6..2dd61e0d3c75714c63cad2f65727692f3430fd2a 100644 --- a/MeshLib/Elements/QuadRule4.cpp +++ b/MeshLib/Elements/QuadRule4.cpp @@ -11,23 +11,21 @@ #include "QuadRule4.h" #include "MathLib/GeometricBasics.h" - #include "MeshLib/Node.h" -namespace MeshLib { - -const unsigned QuadRule4::edge_nodes[4][2] = +namespace MeshLib { - {0, 1}, // Edge 0 - {1, 2}, // Edge 1 - {2, 3}, // Edge 2 - {3, 0} // Edge 3 +const unsigned QuadRule4::edge_nodes[4][2] = { + {0, 1}, // Edge 0 + {1, 2}, // Edge 1 + {2, 3}, // Edge 2 + {3, 0} // Edge 3 }; double QuadRule4::computeVolume(Node const* const* _nodes) { - return MathLib::calcTriangleArea(*_nodes[0], *_nodes[1], *_nodes[2]) - + MathLib::calcTriangleArea(*_nodes[2], *_nodes[3], *_nodes[0]); + return MathLib::calcTriangleArea(*_nodes[0], *_nodes[1], *_nodes[2]) + + MathLib::calcTriangleArea(*_nodes[2], *_nodes[3], *_nodes[0]); } bool QuadRule4::isPntInElement(Node const* const* nodes, @@ -41,7 +39,7 @@ bool QuadRule4::isPntInElement(Node const* const* nodes, unsigned QuadRule4::identifyFace(Node const* const* _nodes, Node* nodes[3]) { - for (unsigned i=0; i<4; i++) + for (unsigned i = 0; i < 4; i++) { unsigned flag(0); for (unsigned j = 0; j < 2; j++) @@ -83,4 +81,4 @@ ElementErrorCode QuadRule4::validate(const Element* e) return error_code; } -} // end namespace MeshLib +} // end namespace MeshLib diff --git a/MeshLib/Elements/QuadRule8.cpp b/MeshLib/Elements/QuadRule8.cpp index 7c4949b8c1c50349bccc9370908cb8317cba93bc..f4c4643bea72252e581a407dbaa0c22c104b3ab1 100644 --- a/MeshLib/Elements/QuadRule8.cpp +++ b/MeshLib/Elements/QuadRule8.cpp @@ -10,14 +10,13 @@ #include "QuadRule8.h" -namespace MeshLib { - -const unsigned QuadRule8::edge_nodes[4][3] = +namespace MeshLib { - {0, 1, 4}, // Edge 0 - {1, 2, 5}, // Edge 1 - {2, 3, 6}, // Edge 2 - {0, 3, 7} // Edge 3 +const unsigned QuadRule8::edge_nodes[4][3] = { + {0, 1, 4}, // Edge 0 + {1, 2, 5}, // Edge 1 + {2, 3, 6}, // Edge 2 + {0, 3, 7} // Edge 3 }; -} // end namespace MeshLib +} // end namespace MeshLib diff --git a/MeshLib/Elements/TemplateElement.cpp b/MeshLib/Elements/TemplateElement.cpp index 5c3237f20b1b9f0118611130df9e679adac996fc..18e766856170b8801e58fc42e783649307ef88cd 100644 --- a/MeshLib/Elements/TemplateElement.cpp +++ b/MeshLib/Elements/TemplateElement.cpp @@ -11,8 +11,8 @@ #include "TemplateElement.h" #include "MeshLib/Elements/Hex.h" -#include "MeshLib/Elements/Point.h" #include "MeshLib/Elements/Line.h" +#include "MeshLib/Elements/Point.h" #include "MeshLib/Elements/Prism.h" #include "MeshLib/Elements/Pyramid.h" #include "MeshLib/Elements/Quad.h" @@ -30,7 +30,7 @@ const unsigned MeshLib::TemplateElement<ELEMENT_RULE>::n_base_nodes; template <class ELEMENT_RULE> const unsigned MeshLib::TemplateElement<ELEMENT_RULE>::dimension; -#endif // WIN32 +#endif // WIN32 template class MeshLib::TemplateElement<MeshLib::HexRule20>; template class MeshLib::TemplateElement<MeshLib::HexRule8>; diff --git a/MeshLib/Elements/TetRule10.cpp b/MeshLib/Elements/TetRule10.cpp index adea61c103f569151c6ab2e03c1e476017c048fc..103929d820ca1857fd6fbb0599b53820afb4800b 100644 --- a/MeshLib/Elements/TetRule10.cpp +++ b/MeshLib/Elements/TetRule10.cpp @@ -13,36 +13,32 @@ #include <array> #include "BaseLib/Logging.h" - #include "MeshLib/Node.h" #include "Tri.h" namespace MeshLib { - -const unsigned TetRule10::face_nodes[4][6] = -{ - {0, 2, 1, 6, 5, 4}, // Face 0 - {0, 1, 3, 4, 8, 7}, // Face 1 - {1, 2, 3, 5, 9, 8}, // Face 2 - {2, 0, 3, 6, 7, 9} // Face 3 +const unsigned TetRule10::face_nodes[4][6] = { + {0, 2, 1, 6, 5, 4}, // Face 0 + {0, 1, 3, 4, 8, 7}, // Face 1 + {1, 2, 3, 5, 9, 8}, // Face 2 + {2, 0, 3, 6, 7, 9} // Face 3 }; -const unsigned TetRule10::edge_nodes[6][3] = -{ - {0, 1, 4}, // Edge 0 - {1, 2, 5}, // Edge 1 - {0, 2, 6}, // Edge 2 - {0, 3, 7}, // Edge 3 - {1, 3, 8}, // Edge 4 - {2, 3, 9} // Edge 5 +const unsigned TetRule10::edge_nodes[6][3] = { + {0, 1, 4}, // Edge 0 + {1, 2, 5}, // Edge 1 + {0, 2, 6}, // Edge 2 + {0, 3, 7}, // Edge 3 + {1, 3, 8}, // Edge 4 + {2, 3, 9} // Edge 5 }; const Element* TetRule10::getFace(const Element* e, unsigned i) { - if (i<n_faces) + if (i < n_faces) { - std::array<Node*,6> nodes{}; + std::array<Node*, 6> nodes{}; for (unsigned j = 0; j < 6; j++) { nodes[j] = const_cast<Node*>(e->getNode(face_nodes[i][j])); @@ -53,4 +49,4 @@ const Element* TetRule10::getFace(const Element* e, unsigned i) return nullptr; } -} // end namespace MeshLib +} // end namespace MeshLib diff --git a/MeshLib/Elements/TetRule4.cpp b/MeshLib/Elements/TetRule4.cpp index 769a26bf72b7cddb862af261904b6e345e6b8901..f43a17ff3859bf24ae6198fc1a65a8f4f3e952dc 100644 --- a/MeshLib/Elements/TetRule4.cpp +++ b/MeshLib/Elements/TetRule4.cpp @@ -13,38 +13,34 @@ #include <array> #include "BaseLib/Logging.h" - +#include "Line.h" #include "MathLib/GeometricBasics.h" - #include "MeshLib/Node.h" #include "Tri.h" -#include "Line.h" - -namespace MeshLib { -const unsigned TetRule4::face_nodes[4][3] = +namespace MeshLib { - {0, 2, 1}, // Face 0 - {0, 1, 3}, // Face 1 - {1, 2, 3}, // Face 2 - {2, 0, 3} // Face 3 +const unsigned TetRule4::face_nodes[4][3] = { + {0, 2, 1}, // Face 0 + {0, 1, 3}, // Face 1 + {1, 2, 3}, // Face 2 + {2, 0, 3} // Face 3 }; -const unsigned TetRule4::edge_nodes[6][2] = -{ - {0, 1}, // Edge 0 - {1, 2}, // Edge 1 - {0, 2}, // Edge 2 - {0, 3}, // Edge 3 - {1, 3}, // Edge 4 - {2, 3} // Edge 5 +const unsigned TetRule4::edge_nodes[6][2] = { + {0, 1}, // Edge 0 + {1, 2}, // Edge 1 + {0, 2}, // Edge 2 + {0, 3}, // Edge 3 + {1, 3}, // Edge 4 + {2, 3} // Edge 5 }; const Element* TetRule4::getFace(const Element* e, unsigned i) { - if (i<n_faces) + if (i < n_faces) { - std::array<Node*,3> nodes{}; + std::array<Node*, 3> nodes{}; for (unsigned j = 0; j < 3; j++) { nodes[j] = const_cast<Node*>(e->getNode(face_nodes[i][j])); @@ -57,7 +53,8 @@ const Element* TetRule4::getFace(const Element* e, unsigned i) double TetRule4::computeVolume(Node const* const* _nodes) { - return MathLib::calcTetrahedronVolume(*_nodes[0], *_nodes[1], *_nodes[2], *_nodes[3]); + return MathLib::calcTetrahedronVolume(*_nodes[0], *_nodes[1], *_nodes[2], + *_nodes[3]); } bool TetRule4::isPntInElement(Node const* const* nodes, @@ -69,7 +66,7 @@ bool TetRule4::isPntInElement(Node const* const* nodes, unsigned TetRule4::identifyFace(Node const* const* _nodes, Node* nodes[3]) { - for (unsigned i=0; i<4; i++) + for (unsigned i = 0; i < 4; i++) { unsigned flag(0); for (unsigned j = 0; j < 3; j++) @@ -94,8 +91,8 @@ ElementErrorCode TetRule4::validate(const Element* e) { ElementErrorCode error_code; error_code[ElementErrorFlag::ZeroVolume] = hasZeroVolume(*e); - error_code[ElementErrorFlag::NodeOrder] = !e->testElementNodeOrder(); + error_code[ElementErrorFlag::NodeOrder] = !e->testElementNodeOrder(); return error_code; } -} // end namespace MeshLib +} // end namespace MeshLib diff --git a/MeshLib/Elements/TriRule3.cpp b/MeshLib/Elements/TriRule3.cpp index 0d363301128f262011c7f7770f67f416cfab86ec..a7fb1c9814f94946214b8131e969eff1e48e8651 100644 --- a/MeshLib/Elements/TriRule3.cpp +++ b/MeshLib/Elements/TriRule3.cpp @@ -11,16 +11,14 @@ #include "TriRule3.h" #include "MathLib/GeometricBasics.h" - #include "MeshLib/Node.h" -namespace MeshLib { - -const unsigned TriRule3::edge_nodes[3][2] = +namespace MeshLib { - {0, 1}, // Edge 0 - {1, 2}, // Edge 1 - {2, 0}, // Edge 2 +const unsigned TriRule3::edge_nodes[3][2] = { + {0, 1}, // Edge 0 + {1, 2}, // Edge 1 + {2, 0}, // Edge 2 }; double TriRule3::computeVolume(Node const* const* _nodes) @@ -37,7 +35,7 @@ bool TriRule3::isPntInElement(Node const* const* nodes, unsigned TriRule3::identifyFace(Node const* const* _nodes, Node* nodes[3]) { - for (unsigned i=0; i<3; i++) + for (unsigned i = 0; i < 3; i++) { unsigned flag(0); for (unsigned j = 0; j < 2; j++) @@ -62,8 +60,8 @@ ElementErrorCode TriRule3::validate(const Element* e) { ElementErrorCode error_code; error_code[ElementErrorFlag::ZeroVolume] = hasZeroVolume(*e); - error_code[ElementErrorFlag::NodeOrder] = !e->testElementNodeOrder(); + error_code[ElementErrorFlag::NodeOrder] = !e->testElementNodeOrder(); return error_code; } -} // end namespace MeshLib +} // end namespace MeshLib diff --git a/MeshLib/Elements/TriRule6.cpp b/MeshLib/Elements/TriRule6.cpp index da31e180f20cdab609ca1f3358725efa5226d070..8c0f1bcd5f4f62eaf05e112d6c8811114e958456 100644 --- a/MeshLib/Elements/TriRule6.cpp +++ b/MeshLib/Elements/TriRule6.cpp @@ -10,13 +10,12 @@ #include "TriRule6.h" -namespace MeshLib { - -const unsigned TriRule6::edge_nodes[3][3] = +namespace MeshLib { - {0, 1, 3}, // Edge 0 - {1, 2, 4}, // Edge 1 - {2, 0, 5}, // Edge 2 +const unsigned TriRule6::edge_nodes[3][3] = { + {0, 1, 3}, // Edge 0 + {1, 2, 4}, // Edge 1 + {2, 0, 5}, // Edge 2 }; -} // end namespace MeshLib +} // end namespace MeshLib diff --git a/MeshLib/IO/Legacy/MeshIO.cpp b/MeshLib/IO/Legacy/MeshIO.cpp index d42080c38ce5f2f0de7c5e2a8e062d96a1ace253..c0b77823e960abae9407c0fca1d5be9965b1be69 100644 --- a/MeshLib/IO/Legacy/MeshIO.cpp +++ b/MeshLib/IO/Legacy/MeshIO.cpp @@ -22,11 +22,9 @@ #include <memory> #include <sstream> -#include "BaseLib/Logging.h" - #include "BaseLib/FileTools.h" +#include "BaseLib/Logging.h" #include "BaseLib/StringTools.h" - #include "MeshLib/Elements/Elements.h" #include "MeshLib/Location.h" #include "MeshLib/Node.h" @@ -36,14 +34,15 @@ namespace MeshLib { namespace IO { -namespace Legacy { +namespace Legacy +{ MeshIO::MeshIO() = default; MeshLib::Mesh* MeshIO::loadMeshFromFile(const std::string& file_name) { INFO("Reading OGS legacy mesh ... "); - std::ifstream in (file_name.c_str(),std::ios::in); + std::ifstream in(file_name.c_str(), std::ios::in); if (!in.is_open()) { WARN("MeshIO::loadMeshFromFile() - Could not open file {:s}.", @@ -58,7 +57,7 @@ MeshLib::Mesh* MeshIO::loadMeshFromFile(const std::string& file_name) std::vector<MeshLib::Element*> elements; std::vector<std::size_t> materials; - if(line_string.find("#FEM_MSH") != std::string::npos) // OGS mesh file + if (line_string.find("#FEM_MSH") != std::string::npos) // OGS mesh file { while (!in.eof()) { @@ -104,17 +103,18 @@ MeshLib::Mesh* MeshIO::loadMeshFromFile(const std::string& file_name) getline(in, line_string); std::stringstream ss(line_string); materials.push_back(readMaterialID(ss)); - MeshLib::Element *elem(readElement(ss,nodes)); - if (elem == nullptr) { + MeshLib::Element* elem(readElement(ss, nodes)); + if (elem == nullptr) + { ERR("Reading mesh element {:d} from file '{:s}' " "failed.", i, file_name); // clean up the elements vector std::for_each(elements.begin(), elements.end(), - std::default_delete<MeshLib::Element>()); + std::default_delete<MeshLib::Element>()); // clean up the nodes vector std::for_each(nodes.begin(), nodes.end(), - std::default_delete<MeshLib::Node>()); + std::default_delete<MeshLib::Node>()); return nullptr; } elements.push_back(elem); @@ -124,7 +124,8 @@ MeshLib::Mesh* MeshIO::loadMeshFromFile(const std::string& file_name) if (elements.empty()) { - ERR ("MeshIO::loadMeshFromFile() - File did not contain element information."); + ERR("MeshIO::loadMeshFromFile() - File did not contain element " + "information."); for (auto& node : nodes) { delete node; @@ -132,8 +133,9 @@ MeshLib::Mesh* MeshIO::loadMeshFromFile(const std::string& file_name) return nullptr; } - MeshLib::Mesh* mesh (new MeshLib::Mesh(BaseLib::extractBaseNameWithoutExtension( - file_name), nodes, elements)); + MeshLib::Mesh* mesh(new MeshLib::Mesh( + BaseLib::extractBaseNameWithoutExtension(file_name), nodes, + elements)); auto* const material_ids = mesh->getProperties().createNewPropertyVector<int>( @@ -159,7 +161,7 @@ MeshLib::Mesh* MeshIO::loadMeshFromFile(const std::string& file_name) return nullptr; } -std::size_t MeshIO::readMaterialID(std::istream & in) +std::size_t MeshIO::readMaterialID(std::istream& in) { unsigned index; unsigned material_id; @@ -170,13 +172,14 @@ std::size_t MeshIO::readMaterialID(std::istream & in) return material_id; } -MeshLib::Element* MeshIO::readElement(std::istream& in, - const std::vector<MeshLib::Node*> &nodes) const +MeshLib::Element* MeshIO::readElement( + std::istream& in, const std::vector<MeshLib::Node*>& nodes) const { std::string elem_type_str; - MeshLib::MeshElemType elem_type (MeshLib::MeshElemType::INVALID); + MeshLib::MeshElemType elem_type(MeshLib::MeshElemType::INVALID); - do { + do + { if (!(in >> elem_type_str)) { return nullptr; @@ -187,134 +190,142 @@ MeshLib::Element* MeshIO::readElement(std::istream& in, auto* idx = new unsigned[8]; MeshLib::Element* elem; - switch(elem_type) + switch (elem_type) { - case MeshLib::MeshElemType::LINE: { - for (int i = 0; i < 2; ++i) + case MeshLib::MeshElemType::LINE: { - if (!(in >> idx[i])) + for (int i = 0; i < 2; ++i) { - return nullptr; + if (!(in >> idx[i])) + { + return nullptr; + } } - } - // edge_nodes array will be deleted from Line object - auto** edge_nodes = new MeshLib::Node*[2]; - for (unsigned k(0); k < 2; ++k) - { - edge_nodes[k] = nodes[idx[k]]; - } - elem = new MeshLib::Line(edge_nodes); - break; - } - case MeshLib::MeshElemType::TRIANGLE: { - for (int i = 0; i < 3; ++i) - { - if (!(in >> idx[i])) + // edge_nodes array will be deleted from Line object + auto** edge_nodes = new MeshLib::Node*[2]; + for (unsigned k(0); k < 2; ++k) { - return nullptr; + edge_nodes[k] = nodes[idx[k]]; } + elem = new MeshLib::Line(edge_nodes); + break; } - auto** tri_nodes = new MeshLib::Node*[3]; - for (unsigned k(0); k < 3; ++k) - { - tri_nodes[k] = nodes[idx[k]]; - } - elem = new MeshLib::Tri(tri_nodes); - break; - } - case MeshLib::MeshElemType::QUAD: { - for (int i = 0; i < 4; ++i) + case MeshLib::MeshElemType::TRIANGLE: { - if (!(in >> idx[i])) + for (int i = 0; i < 3; ++i) { - return nullptr; + if (!(in >> idx[i])) + { + return nullptr; + } } - } - auto** quad_nodes = new MeshLib::Node*[4]; - for (unsigned k(0); k < 4; ++k) - { - quad_nodes[k] = nodes[idx[k]]; - } - elem = new MeshLib::Quad(quad_nodes); - break; - } - case MeshLib::MeshElemType::TETRAHEDRON: { - for (int i = 0; i < 4; ++i) - { - if (!(in >> idx[i])) + auto** tri_nodes = new MeshLib::Node*[3]; + for (unsigned k(0); k < 3; ++k) { - return nullptr; + tri_nodes[k] = nodes[idx[k]]; } + elem = new MeshLib::Tri(tri_nodes); + break; } - auto** tet_nodes = new MeshLib::Node*[4]; - for (unsigned k(0); k < 4; ++k) - { - tet_nodes[k] = nodes[idx[k]]; - } - elem = new MeshLib::Tet(tet_nodes); - break; - } - case MeshLib::MeshElemType::HEXAHEDRON: { - for (int i = 0; i < 8; ++i) + case MeshLib::MeshElemType::QUAD: { - if (!(in >> idx[i])) + for (int i = 0; i < 4; ++i) { - return nullptr; + if (!(in >> idx[i])) + { + return nullptr; + } } + auto** quad_nodes = new MeshLib::Node*[4]; + for (unsigned k(0); k < 4; ++k) + { + quad_nodes[k] = nodes[idx[k]]; + } + elem = new MeshLib::Quad(quad_nodes); + break; } - auto** hex_nodes = new MeshLib::Node*[8]; - for (unsigned k(0); k < 8; ++k) - { - hex_nodes[k] = nodes[idx[k]]; - } - elem = new MeshLib::Hex(hex_nodes); - break; - } - case MeshLib::MeshElemType::PYRAMID: { - for (int i = 0; i < 5; ++i) + case MeshLib::MeshElemType::TETRAHEDRON: { - if (!(in >> idx[i])) + for (int i = 0; i < 4; ++i) + { + if (!(in >> idx[i])) + { + return nullptr; + } + } + auto** tet_nodes = new MeshLib::Node*[4]; + for (unsigned k(0); k < 4; ++k) { - return nullptr; + tet_nodes[k] = nodes[idx[k]]; } + elem = new MeshLib::Tet(tet_nodes); + break; } - auto** pyramid_nodes = new MeshLib::Node*[5]; - for (unsigned k(0); k < 5; ++k) + case MeshLib::MeshElemType::HEXAHEDRON: { - pyramid_nodes[k] = nodes[idx[k]]; + for (int i = 0; i < 8; ++i) + { + if (!(in >> idx[i])) + { + return nullptr; + } + } + auto** hex_nodes = new MeshLib::Node*[8]; + for (unsigned k(0); k < 8; ++k) + { + hex_nodes[k] = nodes[idx[k]]; + } + elem = new MeshLib::Hex(hex_nodes); + break; } - elem = new MeshLib::Pyramid(pyramid_nodes); - break; - } - case MeshLib::MeshElemType::PRISM: { - for (int i = 0; i < 6; ++i) + case MeshLib::MeshElemType::PYRAMID: { - if (!(in >> idx[i])) + for (int i = 0; i < 5; ++i) { - return nullptr; + if (!(in >> idx[i])) + { + return nullptr; + } + } + auto** pyramid_nodes = new MeshLib::Node*[5]; + for (unsigned k(0); k < 5; ++k) + { + pyramid_nodes[k] = nodes[idx[k]]; } + elem = new MeshLib::Pyramid(pyramid_nodes); + break; } - auto** prism_nodes = new MeshLib::Node*[6]; - for (unsigned k(0); k < 6; ++k) + case MeshLib::MeshElemType::PRISM: { - prism_nodes[k] = nodes[idx[k]]; + for (int i = 0; i < 6; ++i) + { + if (!(in >> idx[i])) + { + return nullptr; + } + } + auto** prism_nodes = new MeshLib::Node*[6]; + for (unsigned k(0); k < 6; ++k) + { + prism_nodes[k] = nodes[idx[k]]; + } + elem = new MeshLib::Prism(prism_nodes); + break; } - elem = new MeshLib::Prism(prism_nodes); - break; - } - default: - elem = nullptr; - break; + default: + elem = nullptr; + break; } - delete [] idx; + delete[] idx; return elem; } bool MeshIO::write() { - if(!_mesh) { + if (!_mesh) + { WARN("MeshIO::write(): Cannot write: no mesh object specified."); return false; } @@ -326,7 +337,8 @@ bool MeshIO::write() << " "; const std::size_t n_nodes(_mesh->getNumberOfNodes()); out << n_nodes << "\n"; - for (std::size_t i(0); i < n_nodes; ++i) { + for (std::size_t i(0); i < n_nodes; ++i) + { out << i << " " << *(_mesh->getNode(i)) << "\n"; } @@ -358,10 +370,11 @@ void MeshIO::writeElements( MeshLib::PropertyVector<int> const* const material_ids, std::ostream& out) const { - const std::size_t ele_vector_size (ele_vec.size()); + const std::size_t ele_vector_size(ele_vec.size()); out << ele_vector_size << "\n"; - for (std::size_t i(0); i < ele_vector_size; ++i) { + for (std::size_t i(0); i < ele_vector_size; ++i) + { out << i << " "; if (!material_ids) { @@ -372,7 +385,7 @@ void MeshIO::writeElements( out << (*material_ids)[i] << " "; } out << this->ElemType2StringOutput(ele_vec[i]->getGeomType()) << " "; - unsigned nElemNodes (ele_vec[i]->getNumberOfBaseNodes()); + unsigned nElemNodes(ele_vec[i]->getNumberOfBaseNodes()); for (std::size_t j = 0; j < nElemNodes; ++j) { out << ele_vec[i]->getNode(j)->getID() << " "; @@ -414,6 +427,6 @@ std::string MeshIO::ElemType2StringOutput(const MeshLib::MeshElemType t) return "none"; } -} // end namespace Legacy -} // end namespace IO +} // end namespace Legacy +} // end namespace IO } // namespace MeshLib diff --git a/MeshLib/IO/MPI_IO/NodePartitionedMeshReader.cpp b/MeshLib/IO/MPI_IO/NodePartitionedMeshReader.cpp index da59dafdb821f26aaa914848613ea0f9429a70f2..ce247aa73e54598697d359a9e297436d020dd9c2 100644 --- a/MeshLib/IO/MPI_IO/NodePartitionedMeshReader.cpp +++ b/MeshLib/IO/MPI_IO/NodePartitionedMeshReader.cpp @@ -29,8 +29,7 @@ // Check if the value can by converted to given type without overflow. template <typename VALUE, typename TYPE> -bool -is_safely_convertable(VALUE const& value) +bool is_safely_convertable(VALUE const& value) { bool const result = value <= std::numeric_limits<TYPE>::max(); if (!result) @@ -67,12 +66,13 @@ void NodePartitionedMeshReader::registerNodeDataMpiType() MPI_Datatype types[count] = {MPI_UNSIGNED_LONG, MPI_DOUBLE}; MPI_Aint displacements[count] = {0, sizeof(NodeData::index)}; - MPI_Type_create_struct(count, blocks, displacements, types, &_mpi_node_type); + MPI_Type_create_struct(count, blocks, displacements, types, + &_mpi_node_type); MPI_Type_commit(&_mpi_node_type); } MeshLib::NodePartitionedMesh* NodePartitionedMeshReader::read( - const std::string &file_name_base) + const std::string& file_name_base) { BaseLib::RunTime timer; timer.start(); @@ -81,9 +81,9 @@ MeshLib::NodePartitionedMesh* NodePartitionedMeshReader::read( // Always try binary file first std::string const fname_new = file_name_base + "_partitioned_msh_cfg" + - std::to_string(_mpi_comm_size) + ".bin"; + std::to_string(_mpi_comm_size) + ".bin"; - if(!BaseLib::IsFileExisting(fname_new)) // binary file does not exist. + if (!BaseLib::IsFileExisting(fname_new)) // binary file does not exist. { OGS_FATAL( "Required binary file {:s} does not exist.\n" @@ -103,9 +103,10 @@ MeshLib::NodePartitionedMesh* NodePartitionedMeshReader::read( } template <typename DATA> -bool -NodePartitionedMeshReader::readDataFromFile(std::string const& filename, - MPI_Offset offset, MPI_Datatype type, DATA& data) const +bool NodePartitionedMeshReader::readDataFromFile(std::string const& filename, + MPI_Offset offset, + MPI_Datatype type, + DATA& data) const { // Check container size if (!is_safely_convertable<std::size_t, int>(data.size())) @@ -118,10 +119,10 @@ NodePartitionedMeshReader::readDataFromFile(std::string const& filename, MPI_File file; char* filename_char = const_cast<char*>(filename.data()); - int const file_status = MPI_File_open(_mpi_comm, filename_char, - MPI_MODE_RDONLY, MPI_INFO_NULL, &file); + int const file_status = MPI_File_open( + _mpi_comm, filename_char, MPI_MODE_RDONLY, MPI_INFO_NULL, &file); - if(file_status != 0) + if (file_status != 0) { ERR("Error opening file {:s}. MPI error code {:d}", filename, file_status); @@ -133,24 +134,25 @@ NodePartitionedMeshReader::readDataFromFile(std::string const& filename, MPI_File_set_view(file, offset, type, type, file_mode, MPI_INFO_NULL); // The static cast is checked above. MPI_File_read(file, data.data(), static_cast<int>(data.size()), type, - MPI_STATUS_IGNORE); + MPI_STATUS_IGNORE); MPI_File_close(&file); return true; } MeshLib::NodePartitionedMesh* NodePartitionedMeshReader::readMesh( - const std::string &file_name_base) + const std::string& file_name_base) { //---------------------------------------------------------------------------------- // Read headers - const std::string fname_header = file_name_base + "_partitioned_msh_"; + const std::string fname_header = file_name_base + "_partitioned_msh_"; const std::string fname_num_p_ext = std::to_string(_mpi_comm_size) + ".bin"; - if (!readDataFromFile(fname_header + "cfg" + fname_num_p_ext, - static_cast<MPI_Offset>( - static_cast<unsigned>(_mpi_rank) * sizeof(_mesh_info)), - MPI_LONG, _mesh_info)) + if (!readDataFromFile( + fname_header + "cfg" + fname_num_p_ext, + static_cast<MPI_Offset>(static_cast<unsigned>(_mpi_rank) * + sizeof(_mesh_info)), + MPI_LONG, _mesh_info)) return nullptr; //---------------------------------------------------------------------------------- @@ -158,7 +160,8 @@ MeshLib::NodePartitionedMesh* NodePartitionedMeshReader::readMesh( std::vector<NodeData> nodes(_mesh_info.nodes); if (!readDataFromFile(fname_header + "nod" + fname_num_p_ext, - static_cast<MPI_Offset>(_mesh_info.offset[2]), _mpi_node_type, nodes)) + static_cast<MPI_Offset>(_mesh_info.offset[2]), + _mpi_node_type, nodes)) return nullptr; std::vector<MeshLib::Node*> mesh_nodes; @@ -167,23 +170,25 @@ MeshLib::NodePartitionedMesh* NodePartitionedMeshReader::readMesh( //---------------------------------------------------------------------------------- // Read non-ghost elements - std::vector<unsigned long> elem_data( - _mesh_info.regular_elements + _mesh_info.offset[0]); + std::vector<unsigned long> elem_data(_mesh_info.regular_elements + + _mesh_info.offset[0]); if (!readDataFromFile(fname_header + "ele" + fname_num_p_ext, - static_cast<MPI_Offset>(_mesh_info.offset[3]), MPI_LONG, elem_data)) + static_cast<MPI_Offset>(_mesh_info.offset[3]), + MPI_LONG, elem_data)) return nullptr; - std::vector<MeshLib::Element*> mesh_elems( - _mesh_info.regular_elements + _mesh_info.ghost_elements); + std::vector<MeshLib::Element*> mesh_elems(_mesh_info.regular_elements + + _mesh_info.ghost_elements); setElements(mesh_nodes, elem_data, mesh_elems); //---------------------------------------------------------------------------------- - //Read ghost element - std::vector<unsigned long> ghost_elem_data( - _mesh_info.ghost_elements + _mesh_info.offset[1]); + // Read ghost element + std::vector<unsigned long> ghost_elem_data(_mesh_info.ghost_elements + + _mesh_info.offset[1]); if (!readDataFromFile(fname_header + "ele_g" + fname_num_p_ext, - static_cast<MPI_Offset>(_mesh_info.offset[4]), MPI_LONG, ghost_elem_data)) + static_cast<MPI_Offset>(_mesh_info.offset[4]), + MPI_LONG, ghost_elem_data)) return nullptr; const bool process_ghost = true; @@ -226,7 +231,8 @@ void NodePartitionedMeshReader::readProperties( return; } std::size_t number_of_properties = 0; - is.read(reinterpret_cast<char*>(&number_of_properties), sizeof(std::size_t)); + is.read(reinterpret_cast<char*>(&number_of_properties), + sizeof(std::size_t)); std::vector<std::optional<MeshLib::IO::PropertyVectorMetaData>> vec_pvmd( number_of_properties); for (std::size_t i(0); i < number_of_properties; ++i) @@ -381,14 +387,15 @@ MeshLib::NodePartitionedMesh* NodePartitionedMeshReader::newMesh( _mesh_info.active_nodes); } -void NodePartitionedMeshReader::setNodes(const std::vector<NodeData> &node_data, - std::vector<MeshLib::Node*> &mesh_node, - std::vector<unsigned long> &glb_node_ids) const +void NodePartitionedMeshReader::setNodes( + const std::vector<NodeData>& node_data, + std::vector<MeshLib::Node*>& mesh_node, + std::vector<unsigned long>& glb_node_ids) const { mesh_node.resize(_mesh_info.nodes); glb_node_ids.resize(_mesh_info.nodes); - for(std::size_t i = 0; i < mesh_node.size(); i++) + for (std::size_t i = 0; i < mesh_node.size(); i++) { NodeData const& nd = node_data[i]; glb_node_ids[i] = nd.index; @@ -397,9 +404,9 @@ void NodePartitionedMeshReader::setNodes(const std::vector<NodeData> &node_data, } void NodePartitionedMeshReader::setElements( - const std::vector<MeshLib::Node*> &mesh_nodes, - const std::vector<unsigned long> &elem_data, - std::vector<MeshLib::Element*> &mesh_elems, const bool ghost) const + const std::vector<MeshLib::Node*>& mesh_nodes, + const std::vector<unsigned long>& elem_data, + std::vector<MeshLib::Element*>& mesh_elems, const bool ghost) const { // Number of elements, either ghost or regular unsigned long const ne = @@ -407,7 +414,7 @@ void NodePartitionedMeshReader::setElements( unsigned long const id_offset_ghost = ghost ? _mesh_info.regular_elements : 0; - for(unsigned long i = 0; i < ne; i++) + for (unsigned long i = 0; i < ne; i++) { unsigned long id_offset_elem = elem_data[i]; @@ -421,8 +428,8 @@ void NodePartitionedMeshReader::setElements( unsigned long const nnodes = elem_data[id_offset_elem++]; MeshLib::Node** elem_nodes = new MeshLib::Node*[nnodes]; - for(unsigned long k = 0; k < nnodes; k++) - elem_nodes[k] = mesh_nodes[ elem_data[id_offset_elem++] ]; + for (unsigned long k = 0; k < nnodes; k++) + elem_nodes[k] = mesh_nodes[elem_data[id_offset_elem++]]; // The element types below are defined by the MeshLib::CellType. switch (static_cast<CellType>(e_type)) @@ -505,5 +512,5 @@ void NodePartitionedMeshReader::setElements( } } } -} // namespace IO -} // namespace MeshLib +} // namespace IO +} // namespace MeshLib diff --git a/MeshLib/IO/VtkIO/PVDFile.cpp b/MeshLib/IO/VtkIO/PVDFile.cpp index 5ba3868dd53c7362de4136644adcf678eceb8ccd..f2be2d4c6bc6d4e5f0088433ef90efdab41ff5e2 100644 --- a/MeshLib/IO/VtkIO/PVDFile.cpp +++ b/MeshLib/IO/VtkIO/PVDFile.cpp @@ -13,6 +13,7 @@ #include <fstream> #include <iomanip> #include <limits> + #include "BaseLib/Error.h" #include "MeshLib/IO/VtkIO/VtuInterface.h" diff --git a/MeshLib/IO/XDMF/HdfData.cpp b/MeshLib/IO/XDMF/HdfData.cpp index 86f6809f8ae2a5dfa0574e870e1f5ff0cc7cad82..0103c67a18e1183112103485998fbf549116bf3c 100644 --- a/MeshLib/IO/XDMF/HdfData.cpp +++ b/MeshLib/IO/XDMF/HdfData.cpp @@ -48,8 +48,10 @@ HdfData::HdfData(void const* data_start, std::size_t const size_partitioned_dim, name(name) { auto const& partition_info = getPartitionInfo(size_partitioned_dim); - DBUG("HdfData: The partition of dataset {:s} has dimension {:d} and offset {:d}.", - name, size_partitioned_dim, partition_info.first); + DBUG( + "HdfData: The partition of dataset {:s} has dimension {:d} and offset " + "{:d}.", + name, size_partitioned_dim, partition_info.first); auto const& offset_partitioned_dim = partition_info.first; offsets = {offset_partitioned_dim, 0}; file_space = {partition_info.second, size_tuple}; diff --git a/MeshLib/IO/readMeshFromFile.cpp b/MeshLib/IO/readMeshFromFile.cpp index d0e67283b2515be0feb9c43f3c4700f2878a60fb..6c60dd6dc98dd7a6d35cd1fbcecab73b901c2807 100644 --- a/MeshLib/IO/readMeshFromFile.cpp +++ b/MeshLib/IO/readMeshFromFile.cpp @@ -23,15 +23,12 @@ #include <boost/algorithm/string/erase.hpp> -#include "BaseLib/Logging.h" - #include "BaseLib/FileTools.h" +#include "BaseLib/Logging.h" #include "BaseLib/StringTools.h" - -#include "MeshLib/Mesh.h" - #include "MeshLib/IO/Legacy/MeshIO.h" #include "MeshLib/IO/VtkIO/VtuInterface.h" +#include "MeshLib/Mesh.h" #ifdef USE_PETSC #include "MeshLib/IO/MPI_IO/NodePartitionedMeshReader.h" @@ -42,7 +39,7 @@ namespace MeshLib { namespace IO { -MeshLib::Mesh* readMeshFromFile(const std::string &file_name) +MeshLib::Mesh* readMeshFromFile(const std::string& file_name) { #ifdef USE_PETSC int world_size; @@ -50,14 +47,15 @@ MeshLib::Mesh* readMeshFromFile(const std::string &file_name) if (world_size > 1) { MeshLib::IO::NodePartitionedMeshReader read_pmesh(PETSC_COMM_WORLD); - const std::string file_name_base = BaseLib::dropFileExtension(file_name); + const std::string file_name_base = + BaseLib::dropFileExtension(file_name); return read_pmesh.read(file_name_base); } else if (world_size == 1) { MeshLib::Mesh* mesh = readMeshFromFileSerial(file_name); - MeshLib::NodePartitionedMesh* part_mesh - = new MeshLib::NodePartitionedMesh(*mesh); + MeshLib::NodePartitionedMesh* part_mesh = + new MeshLib::NodePartitionedMesh(*mesh); delete mesh; return part_mesh; } @@ -67,7 +65,7 @@ MeshLib::Mesh* readMeshFromFile(const std::string &file_name) #endif } -MeshLib::Mesh* readMeshFromFileSerial(const std::string &file_name) +MeshLib::Mesh* readMeshFromFileSerial(const std::string& file_name) { if (BaseLib::hasFileExtension(".msh", file_name)) { @@ -90,6 +88,5 @@ MeshLib::Mesh* readMeshFromFileSerial(const std::string &file_name) return nullptr; } - -} // end namespace IO -} // end namespace MeshLib +} // end namespace IO +} // end namespace MeshLib diff --git a/MeshLib/IO/writeMeshToFile.cpp b/MeshLib/IO/writeMeshToFile.cpp index 1adb694ad0d2418fc0d3bdc795d171fe041c3cec..3f899d0d015575f34d7c4f15bb6953ac13e737bc 100644 --- a/MeshLib/IO/writeMeshToFile.cpp +++ b/MeshLib/IO/writeMeshToFile.cpp @@ -44,9 +44,8 @@ int writeMeshToFile(const MeshLib::Mesh& mesh, #ifdef OGS_USE_XDMF if (file_path.extension().string() == ".xdmf") { - auto writer = - std::make_unique<MeshLib::IO::XdmfHdfWriter>(MeshLib::IO::XdmfHdfWriter( - mesh, file_path, 0)); + auto writer = std::make_unique<MeshLib::IO::XdmfHdfWriter>( + MeshLib::IO::XdmfHdfWriter(mesh, file_path, 0)); // \TODO Errorhandling, Change data model into static and time depended, // then is is not necessary to give time step 0 a special treatment // here diff --git a/MeshLib/Location.cpp b/MeshLib/Location.cpp index e8d1c459be58bc5a864127785d5b5fedce857fec..b0617f62a7e3d5c3147fa5bac9ab69ff3fb761bf 100644 --- a/MeshLib/Location.cpp +++ b/MeshLib/Location.cpp @@ -9,30 +9,33 @@ */ #include "Location.h" + #include <ostream> namespace MeshLib { - std::ostream& operator<<(std::ostream& os, MeshItemType const& t) { switch (t) { - case MeshItemType::Node: return os << "N"; - case MeshItemType::Edge: return os << "E"; - case MeshItemType::Face: return os << "F"; - case MeshItemType::Cell: return os << "C"; - case MeshItemType::IntegrationPoint: return os << "I"; + case MeshItemType::Node: + return os << "N"; + case MeshItemType::Edge: + return os << "E"; + case MeshItemType::Face: + return os << "F"; + case MeshItemType::Cell: + return os << "C"; + case MeshItemType::IntegrationPoint: + return os << "I"; }; return os; } std::ostream& operator<<(std::ostream& os, Location const& l) { - return os << "(" << l.mesh_id - << ", " << l.item_type - << ", " << l.item_id - << ")"; + return os << "(" << l.mesh_id << ", " << l.item_type << ", " << l.item_id + << ")"; } -} // namespace MeshLib +} // namespace MeshLib diff --git a/MeshLib/Mesh.cpp b/MeshLib/Mesh.cpp index 35805ad4a6ce41dafe586fbb7f91d5e499e69545..ee5c3ac50c9484084166e45c03010766a128146b 100644 --- a/MeshLib/Mesh.cpp +++ b/MeshLib/Mesh.cpp @@ -19,14 +19,13 @@ #include <utility> #include "BaseLib/RunTime.h" - #include "Elements/Element.h" -#include "Elements/Tri.h" -#include "Elements/Quad.h" -#include "Elements/Tet.h" #include "Elements/Hex.h" -#include "Elements/Pyramid.h" #include "Elements/Prism.h" +#include "Elements/Pyramid.h" +#include "Elements/Quad.h" +#include "Elements/Tet.h" +#include "Elements/Tri.h" /// Mesh counter used to uniquely identify meshes by id. static std::size_t global_mesh_counter = 0; @@ -80,16 +79,16 @@ Mesh::Mesh(const Mesh& mesh) _n_base_nodes(mesh.getNumberOfBaseNodes()), _properties(mesh._properties) { - const std::vector<Node*>& nodes (mesh.getNodes()); - const std::size_t nNodes (nodes.size()); + const std::vector<Node*>& nodes(mesh.getNodes()); + const std::size_t nNodes(nodes.size()); for (unsigned i = 0; i < nNodes; ++i) { _nodes[i] = new Node(*nodes[i]); } - const std::vector<Element*>& elements (mesh.getElements()); - const std::size_t nElements (elements.size()); - for (unsigned i=0; i<nElements; ++i) + const std::vector<Element*>& elements(mesh.getElements()); + const std::size_t nElements(elements.size()); + for (unsigned i = 0; i < nElements; ++i) { const std::size_t nElemNodes = elements[i]->getNumberOfNodes(); _elements[i] = elements[i]->clone(); @@ -104,19 +103,19 @@ Mesh::Mesh(const Mesh& mesh) this->setDimension(); } this->setElementsConnectedToNodes(); - //this->setNodesConnectedByElements(); + // this->setNodesConnectedByElements(); this->setElementNeighbors(); } Mesh::~Mesh() { - const std::size_t nElements (_elements.size()); + const std::size_t nElements(_elements.size()); for (std::size_t i = 0; i < nElements; ++i) { delete _elements[i]; } - const std::size_t nNodes (_nodes.size()); + const std::size_t nNodes(_nodes.size()); for (std::size_t i = 0; i < nNodes; ++i) { delete _nodes[i]; @@ -128,7 +127,7 @@ void Mesh::addElement(Element* elem) _elements.push_back(elem); // add element information to nodes - unsigned nNodes (elem->getNumberOfNodes()); + unsigned nNodes(elem->getNumberOfNodes()); for (unsigned i = 0; i < nNodes; ++i) { elem->_nodes[i]->addElement(elem); @@ -159,7 +158,7 @@ void Mesh::recalculateMaxBaseNodeId() void Mesh::resetElementIDs() { - const std::size_t nElements (this->_elements.size()); + const std::size_t nElements(this->_elements.size()); for (unsigned i = 0; i < nElements; ++i) { _elements[i]->setID(i); @@ -168,7 +167,7 @@ void Mesh::resetElementIDs() void Mesh::setDimension() { - const std::size_t nElements (_elements.size()); + const std::size_t nElements(_elements.size()); for (unsigned i = 0; i < nElements; ++i) { if (_elements[i]->getDimension() > _mesh_dimension) @@ -209,17 +208,22 @@ void Mesh::setElementNeighbors() std::vector<Element*> neighbors; for (auto element : _elements) { - // create vector with all elements connected to current element (includes lots of doubles!) - const std::size_t nNodes (element->getNumberOfBaseNodes()); - for (unsigned n(0); n<nNodes; ++n) + // create vector with all elements connected to current element + // (includes lots of doubles!) + const std::size_t nNodes(element->getNumberOfBaseNodes()); + for (unsigned n(0); n < nNodes; ++n) { - std::vector<Element*> const& conn_elems ((element->getNode(n)->getElements())); - neighbors.insert(neighbors.end(), conn_elems.begin(), conn_elems.end()); + std::vector<Element*> const& conn_elems( + (element->getNode(n)->getElements())); + neighbors.insert(neighbors.end(), conn_elems.begin(), + conn_elems.end()); } std::sort(neighbors.begin(), neighbors.end()); - auto const neighbors_new_end = std::unique(neighbors.begin(), neighbors.end()); + auto const neighbors_new_end = + std::unique(neighbors.begin(), neighbors.end()); - for (auto neighbor = neighbors.begin(); neighbor != neighbors_new_end; ++neighbor) + for (auto neighbor = neighbors.begin(); neighbor != neighbors_new_end; + ++neighbor) { std::optional<unsigned> const opposite_face_id = element->addNeighbor(*neighbor); @@ -257,8 +261,9 @@ void Mesh::setNodesConnectedByElements() // Make nodes unique and sorted by their ids. // This relies on the node's id being equivalent to it's address. std::sort(adjacent_nodes.begin(), adjacent_nodes.end(), - [](Node* a, Node* b) { return a->getID() < b->getID(); }); - auto const last = std::unique(adjacent_nodes.begin(), adjacent_nodes.end()); + [](Node* a, Node* b) { return a->getID() < b->getID(); }); + auto const last = + std::unique(adjacent_nodes.begin(), adjacent_nodes.end()); adjacent_nodes.erase(last, adjacent_nodes.end()); node->setConnectedNodes(adjacent_nodes); @@ -269,7 +274,8 @@ void Mesh::checkNonlinearNodeIDs() const { for (MeshLib::Element const* e : _elements) { - for (unsigned i=e->getNumberOfBaseNodes(); i<e->getNumberOfNodes(); i++) + for (unsigned i = e->getNumberOfBaseNodes(); i < e->getNumberOfNodes(); + i++) { if (e->getNodeIndex(i) >= getNumberOfBaseNodes()) { @@ -288,13 +294,13 @@ void Mesh::checkNonlinearNodeIDs() const bool Mesh::hasNonlinearElement() const { - return std::any_of(std::begin(_elements), std::end(_elements), - [](Element const* const e) { + return std::any_of( + std::begin(_elements), std::end(_elements), [](Element const* const e) { return e->getNumberOfNodes() != e->getNumberOfBaseNodes(); }); } -void scaleMeshPropertyVector(MeshLib::Mesh & mesh, +void scaleMeshPropertyVector(MeshLib::Mesh& mesh, std::string const& property_name, double factor) { @@ -303,7 +309,7 @@ void scaleMeshPropertyVector(MeshLib::Mesh & mesh, WARN("Did not find PropertyVector '{:s}' for scaling.", property_name); return; } - auto & pv = *mesh.getProperties().getPropertyVector<double>(property_name); + auto& pv = *mesh.getProperties().getPropertyVector<double>(property_name); std::transform(pv.begin(), pv.end(), pv.begin(), [factor](auto const& v) { return v * factor; }); } diff --git a/MeshLib/MeshEditing/AddLayerToMesh.cpp b/MeshLib/MeshEditing/AddLayerToMesh.cpp index 89af403f11cbe72292416582515a62557d28e299..5ad8c4e2c098106b80093fa978e29e2e64d8f8ad 100644 --- a/MeshLib/MeshEditing/AddLayerToMesh.cpp +++ b/MeshLib/MeshEditing/AddLayerToMesh.cpp @@ -19,17 +19,15 @@ #include <vector> #include "BaseLib/Logging.h" - -#include "MeshLib/Mesh.h" -#include "MeshLib/Node.h" #include "MeshLib/Elements/Elements.h" -#include "MeshLib/MeshSurfaceExtraction.h" +#include "MeshLib/Mesh.h" #include "MeshLib/MeshEditing/DuplicateMeshComponents.h" #include "MeshLib/MeshEditing/FlipElements.h" +#include "MeshLib/MeshSurfaceExtraction.h" +#include "MeshLib/Node.h" namespace MeshLib { - /** Extrudes point, line, triangle or quad elements to its higher dimensional * versions, i.e. line, quad, prism, hexahedron. * @@ -42,8 +40,9 @@ namespace MeshLib * * @return extruded element (point -> line, line -> quad, tri -> prism, quad -> * hexahedron) -*/ -MeshLib::Element* extrudeElement(std::vector<MeshLib::Node*> const& subsfc_nodes, + */ +MeshLib::Element* extrudeElement( + std::vector<MeshLib::Node*> const& subsfc_nodes, MeshLib::Element const& sfc_elem, MeshLib::PropertyVector<std::size_t> const& sfc_to_subsfc_id_map, std::map<std::size_t, std::size_t> const& subsfc_sfc_id_map) @@ -58,12 +57,12 @@ MeshLib::Element* extrudeElement(std::vector<MeshLib::Node*> const& subsfc_nodes new MeshLib::Node*[2 * nElemNodes] }; - for (unsigned j=0; j<nElemNodes; ++j) + for (unsigned j = 0; j < nElemNodes; ++j) { std::size_t const subsfc_id( sfc_to_subsfc_id_map[sfc_elem.getNode(j)->getID()]); new_nodes[j] = subsfc_nodes[subsfc_id]; - std::size_t new_idx = (nElemNodes==2) ? (3-j) : (nElemNodes+j); + std::size_t new_idx = (nElemNodes == 2) ? (3 - j) : (nElemNodes + j); new_nodes[new_idx] = subsfc_nodes[subsfc_sfc_id_map.at(subsfc_id)]; } @@ -91,7 +90,7 @@ MeshLib::Mesh* addLayerToMesh(MeshLib::Mesh const& mesh, double thickness, double const flag = on_top ? -1.0 : 1.0; Eigen::Vector3d const dir({0, 0, flag}); double const angle(90); - std::unique_ptr<MeshLib::Mesh> sfc_mesh (nullptr); + std::unique_ptr<MeshLib::Mesh> sfc_mesh(nullptr); std::string const prop_name("bulk_node_ids"); @@ -210,4 +209,4 @@ MeshLib::Mesh* addLayerToMesh(MeshLib::Mesh const& mesh, double thickness, return new_mesh; } -} // namespace MeshLib +} // namespace MeshLib diff --git a/MeshLib/MeshEditing/DuplicateMeshComponents.cpp b/MeshLib/MeshEditing/DuplicateMeshComponents.cpp index fc4863a5c29dffc2f24c89ea6b56c1d0a880bb5f..14b4884c450c77aab57171f22a26a8990beb072c 100644 --- a/MeshLib/MeshEditing/DuplicateMeshComponents.cpp +++ b/MeshLib/MeshEditing/DuplicateMeshComponents.cpp @@ -14,8 +14,8 @@ #include "DuplicateMeshComponents.h" -#include "MeshLib/Mesh.h" #include "MeshLib/Elements/Elements.h" +#include "MeshLib/Mesh.h" namespace MeshLib { @@ -74,10 +74,9 @@ MeshLib::Element* copyElement(MeshLib::Element const* const element, return new E(new_nodes); } -MeshLib::Element* copyElement( - MeshLib::Element const* const element, - const std::vector<MeshLib::Node*>& nodes, - std::vector<std::size_t> const* const id_map) +MeshLib::Element* copyElement(MeshLib::Element const* const element, + const std::vector<MeshLib::Node*>& nodes, + std::vector<std::size_t> const* const id_map) { switch (element->getCellType()) { diff --git a/MeshLib/MeshEditing/ElementValueModification.cpp b/MeshLib/MeshEditing/ElementValueModification.cpp index 4a727c233eade11a1f0f44d531b8f35666cbf4b3..bbebbc9c7323ae0d6816b80b8a16fb48cd968b46 100644 --- a/MeshLib/MeshEditing/ElementValueModification.cpp +++ b/MeshLib/MeshEditing/ElementValueModification.cpp @@ -17,16 +17,16 @@ #include <algorithm> #include "BaseLib/Logging.h" - #include "MeshLib/Elements/Element.h" #include "MeshLib/Mesh.h" #include "MeshLib/PropertyVector.h" -namespace MeshLib { - -bool ElementValueModification::replace(MeshLib::Mesh &mesh, - std::string const& property_name, int const old_value, int const new_value, - bool replace_if_exists) +namespace MeshLib +{ +bool ElementValueModification::replace(MeshLib::Mesh& mesh, + std::string const& property_name, + int const old_value, int const new_value, + bool replace_if_exists) { MeshLib::PropertyVector<int>* property_value_vec = nullptr; try @@ -70,13 +70,15 @@ bool ElementValueModification::replace(MeshLib::Mesh &mesh, return true; } -bool ElementValueModification::replace(MeshLib::Mesh &mesh, - int const old_value, int const new_value, bool replace_if_exists) +bool ElementValueModification::replace(MeshLib::Mesh& mesh, int const old_value, + int const new_value, + bool replace_if_exists) { - return replace(mesh, "MaterialIDs", old_value, new_value, replace_if_exists); + return replace(mesh, "MaterialIDs", old_value, new_value, + replace_if_exists); } -std::size_t ElementValueModification::condense(MeshLib::Mesh &mesh) +std::size_t ElementValueModification::condense(MeshLib::Mesh& mesh) { MeshLib::PropertyVector<int>* property_value_vector = nullptr; try @@ -110,7 +112,9 @@ std::size_t ElementValueModification::condense(MeshLib::Mesh &mesh) return nValues; } -std::size_t ElementValueModification::setByElementType(MeshLib::Mesh &mesh, MeshElemType ele_type, int const new_value) +std::size_t ElementValueModification::setByElementType(MeshLib::Mesh& mesh, + MeshElemType ele_type, + int const new_value) { MeshLib::PropertyVector<int>* property_value_vector = nullptr; try @@ -139,4 +143,4 @@ std::size_t ElementValueModification::setByElementType(MeshLib::Mesh &mesh, Mesh return cnt; } -} // end namespace MeshLib +} // end namespace MeshLib diff --git a/MeshLib/MeshEditing/FlipElements.cpp b/MeshLib/MeshEditing/FlipElements.cpp index b2b3ec1f1f299a41547a51bab370ab47582c7a6b..59cf46ee38ed8fefc19ff423f22be0034a98d72d 100644 --- a/MeshLib/MeshEditing/FlipElements.cpp +++ b/MeshLib/MeshEditing/FlipElements.cpp @@ -10,12 +10,12 @@ #include "FlipElements.h" +#include "DuplicateMeshComponents.h" #include "MeshLib/Elements/Element.h" -#include "MeshLib/Node.h" #include "MeshLib/Elements/Line.h" -#include "MeshLib/Elements/Tri.h" #include "MeshLib/Elements/Quad.h" -#include "DuplicateMeshComponents.h" +#include "MeshLib/Elements/Tri.h" +#include "MeshLib/Node.h" namespace MeshLib { diff --git a/MeshLib/MeshEditing/Mesh2MeshPropertyInterpolation.cpp b/MeshLib/MeshEditing/Mesh2MeshPropertyInterpolation.cpp index 9729f8da8b24eb9d16c874c4a424041e96af0b65..b8ede148ebc3544cf31b74f5edcb439eae294f32 100644 --- a/MeshLib/MeshEditing/Mesh2MeshPropertyInterpolation.cpp +++ b/MeshLib/MeshEditing/Mesh2MeshPropertyInterpolation.cpp @@ -12,31 +12,31 @@ * */ -#include <vector> -#include <fstream> - #include "Mesh2MeshPropertyInterpolation.h" -#include "BaseLib/Logging.h" +#include <fstream> +#include <vector> +#include "BaseLib/Logging.h" #include "GeoLib/AABB.h" #include "GeoLib/Grid.h" - -#include "MeshLib/MeshEnums.h" +#include "MeshLib/Elements/Element.h" #include "MeshLib/Mesh.h" +#include "MeshLib/MeshEnums.h" #include "MeshLib/Node.h" -#include "MeshLib/Elements/Element.h" - -namespace MeshLib { +namespace MeshLib +{ Mesh2MeshPropertyInterpolation::Mesh2MeshPropertyInterpolation( Mesh const& src_mesh, std::string const& property_name) : _src_mesh(src_mesh), _property_name(property_name) -{} +{ +} bool Mesh2MeshPropertyInterpolation::setPropertiesForMesh(Mesh& dest_mesh) const { - if (_src_mesh.getDimension() != dest_mesh.getDimension()) { + if (_src_mesh.getDimension() != dest_mesh.getDimension()) + { ERR("MeshLib::Mesh2MeshPropertyInterpolation::setPropertiesForMesh() " "dimension of source (dim = {:d}) and destination (dim = {:d}) " "mesh does not match.", @@ -44,7 +44,8 @@ bool Mesh2MeshPropertyInterpolation::setPropertiesForMesh(Mesh& dest_mesh) const return false; } - if (_src_mesh.getDimension() != 2) { + if (_src_mesh.getDimension() != 2) + { WARN( "MeshLib::Mesh2MeshPropertyInterpolation::setPropertiesForMesh() " "implemented only for 2D case at the moment."); @@ -147,8 +148,9 @@ void Mesh2MeshPropertyInterpolation::interpolatePropertiesForMesh( } } -void Mesh2MeshPropertyInterpolation::interpolateElementPropertiesToNodeProperties( - std::vector<double> &interpolated_properties) const +void Mesh2MeshPropertyInterpolation:: + interpolateElementPropertiesToNodeProperties( + std::vector<double>& interpolated_properties) const { // fetch the source of property values if (!_src_mesh.getProperties().existsPropertyVector<double>(_property_name)) @@ -175,4 +177,4 @@ void Mesh2MeshPropertyInterpolation::interpolateElementPropertiesToNodePropertie } } -} // end namespace MeshLib +} // end namespace MeshLib diff --git a/MeshLib/MeshEditing/ProjectPointOnMesh.cpp b/MeshLib/MeshEditing/ProjectPointOnMesh.cpp index 0dce76e23c9beadff2aa7615bd8919d907d52383..88930006e5858debe69454c5bdd90ce73414dcf0 100644 --- a/MeshLib/MeshEditing/ProjectPointOnMesh.cpp +++ b/MeshLib/MeshEditing/ProjectPointOnMesh.cpp @@ -10,7 +10,6 @@ #include "GeoLib/AnalyticalGeometry.h" #include "MathLib/Point3d.h" - #include "MeshLib/Elements/Quad.h" #include "MeshLib/Elements/Tri.h" #include "MeshLib/MeshSearch/MeshElementGrid.h" @@ -19,7 +18,6 @@ namespace MeshLib { namespace ProjectPointOnMesh { - MeshLib::Element const* getProjectedElement( std::vector<const MeshLib::Element*> const& elements, MeshLib::Node const& node) diff --git a/MeshLib/MeshEditing/RasterDataToMesh.cpp b/MeshLib/MeshEditing/RasterDataToMesh.cpp index c2f426286b50a028c843238bae53c4de8f6b44b1..64d7977f3f1e30a2377994a458561a00a5b06a36 100644 --- a/MeshLib/MeshEditing/RasterDataToMesh.cpp +++ b/MeshLib/MeshEditing/RasterDataToMesh.cpp @@ -11,12 +11,11 @@ #include "RasterDataToMesh.h" #include "BaseLib/StringTools.h" -#include "MeshLib/Node.h" #include "MeshLib/Elements/Element.h" +#include "MeshLib/Node.h" namespace MeshLib { - namespace RasterDataToMesh { static bool checkMesh(MeshLib::Mesh const& mesh) diff --git a/MeshLib/MeshEditing/RemoveMeshComponents.cpp b/MeshLib/MeshEditing/RemoveMeshComponents.cpp index e05da29001a7b6980cdef92e71d90603cdd00648..f9fe696ca2b0cbc7ab1f1676cfd520ea409c156a 100644 --- a/MeshLib/MeshEditing/RemoveMeshComponents.cpp +++ b/MeshLib/MeshEditing/RemoveMeshComponents.cpp @@ -10,11 +10,11 @@ #include "RemoveMeshComponents.h" +#include "DuplicateMeshComponents.h" #include "MeshLib/Elements/Element.h" -#include "MeshLib/Node.h" #include "MeshLib/MeshSearch/ElementSearch.h" #include "MeshLib/MeshSearch/NodeSearch.h" -#include "DuplicateMeshComponents.h" +#include "MeshLib/Node.h" namespace MeshLib { @@ -174,4 +174,3 @@ MeshLib::Mesh* removeNodes(const MeshLib::Mesh& mesh, return nullptr; } } // end namespace MeshLib - diff --git a/MeshLib/MeshEnums.cpp b/MeshLib/MeshEnums.cpp index b6dcd5729bf9f7095e6946af2ada6f77c2032481..d36438b0b2bcabe58fb77c09a533d2ffaf2cdbbc 100644 --- a/MeshLib/MeshEnums.cpp +++ b/MeshLib/MeshEnums.cpp @@ -16,7 +16,8 @@ #include <boost/algorithm/string/predicate.hpp> -namespace MeshLib { +namespace MeshLib +{ std::string MeshElemType2String(const MeshElemType t) { if (t == MeshElemType::POINT) @@ -91,7 +92,7 @@ std::string MeshElemType2StringShort(const MeshElemType t) return "none"; } -MeshElemType String2MeshElemType(const std::string &s) +MeshElemType String2MeshElemType(const std::string& s) { if (boost::iequals(s, "point")) { @@ -155,8 +156,8 @@ std::vector<std::string> getMeshElemTypeStringsShort() std::string CellType2String(const CellType t) { -#define RETURN_CELL_TYPE_STR(t, type)\ - if ((t) == CellType::type)\ +#define RETURN_CELL_TYPE_STR(t, type) \ + if ((t) == CellType::type) \ return #type; RETURN_CELL_TYPE_STR(t, POINT1); diff --git a/MeshLib/MeshGenerators/LayeredMeshGenerator.cpp b/MeshLib/MeshGenerators/LayeredMeshGenerator.cpp index c47bbce05c6d68c937d4ee9b30fbcd30110e973c..584d1be408bd35ceb6a42b2e3845b62797ceaa53 100644 --- a/MeshLib/MeshGenerators/LayeredMeshGenerator.cpp +++ b/MeshLib/MeshGenerators/LayeredMeshGenerator.cpp @@ -14,21 +14,19 @@ #include "LayeredMeshGenerator.h" -#include <vector> #include <fstream> +#include <vector> #include "BaseLib/Logging.h" - #include "GeoLib/Raster.h" - -#include "MeshLib/Mesh.h" -#include "MeshLib/Node.h" #include "MeshLib/Elements/Element.h" +#include "MeshLib/Mesh.h" +#include "MeshLib/MeshEditing/RemoveMeshComponents.h" #include "MeshLib/MeshInformation.h" -#include "MeshLib/PropertyVector.h" -#include "MeshLib/Properties.h" #include "MeshLib/MeshSearch/NodeSearch.h" -#include "MeshLib/MeshEditing/RemoveMeshComponents.h" +#include "MeshLib/Node.h" +#include "MeshLib/Properties.h" +#include "MeshLib/PropertyVector.h" bool LayeredMeshGenerator::createLayers( MeshLib::Mesh const& mesh, @@ -41,7 +39,8 @@ bool LayeredMeshGenerator::createLayers( return false; } - auto const elem_count = MeshLib::MeshInformation::getNumberOfElementTypes(mesh); + auto const elem_count = + MeshLib::MeshInformation::getNumberOfElementTypes(mesh); if (elem_count.find(MeshLib::MeshElemType::QUAD) != elem_count.end()) { ERR("Input mesh contains QUAD-elements. Please use input mesh " @@ -49,13 +48,16 @@ bool LayeredMeshGenerator::createLayers( return false; } - bool result = createRasterLayers(mesh, rasters, minimum_thickness, noDataReplacementValue); - std::for_each(rasters.begin(), rasters.end(), [](GeoLib::Raster const*const raster){ delete raster; }); + bool result = createRasterLayers( + mesh, rasters, minimum_thickness, noDataReplacementValue); + std::for_each(rasters.begin(), + rasters.end(), + [](GeoLib::Raster const* const raster) { delete raster; }); return result; } -std::unique_ptr<MeshLib::Mesh> -LayeredMeshGenerator::getMesh(std::string const& mesh_name) const +std::unique_ptr<MeshLib::Mesh> LayeredMeshGenerator::getMesh( + std::string const& mesh_name) const { if (_nodes.empty()) { @@ -86,9 +88,11 @@ LayeredMeshGenerator::getMesh(std::string const& mesh_name) const "element number"); } - std::unique_ptr<MeshLib::Mesh> result(new MeshLib::Mesh(mesh_name, _nodes, _elements, properties)); + std::unique_ptr<MeshLib::Mesh> result( + new MeshLib::Mesh(mesh_name, _nodes, _elements, properties)); MeshLib::NodeSearch ns(*result); - if (ns.searchUnused() > 0) { + if (ns.searchUnused() > 0) + { std::unique_ptr<MeshLib::Mesh> new_mesh( MeshLib::removeNodes(*result, ns.getSearchedNodeIDs(), mesh_name)); return new_mesh; @@ -96,19 +100,22 @@ LayeredMeshGenerator::getMesh(std::string const& mesh_name) const return result; } -double LayeredMeshGenerator::calcEpsilon(GeoLib::Raster const& low, GeoLib::Raster const& high) +double LayeredMeshGenerator::calcEpsilon(GeoLib::Raster const& low, + GeoLib::Raster const& high) { - const double max (*std::max_element(high.begin(), high.end())); - const double min (*std::min_element( low.begin(), low.end())); - return ((max-min)*1e-07); + const double max(*std::max_element(high.begin(), high.end())); + const double min(*std::min_element(low.begin(), low.end())); + return ((max - min) * 1e-07); } -MeshLib::Node* LayeredMeshGenerator::getNewLayerNode(MeshLib::Node const& dem_node, - MeshLib::Node const& last_layer_node, - GeoLib::Raster const& raster, - std::size_t new_node_id) const +MeshLib::Node* LayeredMeshGenerator::getNewLayerNode( + MeshLib::Node const& dem_node, + MeshLib::Node const& last_layer_node, + GeoLib::Raster const& raster, + std::size_t new_node_id) const { - double const elevation = std::min(raster.interpolateValueAtPoint(dem_node), dem_node[2]); + double const elevation = + std::min(raster.interpolateValueAtPoint(dem_node), dem_node[2]); if ((std::abs(elevation - raster.getHeader().no_data) < std::numeric_limits<double>::epsilon()) || diff --git a/MeshLib/MeshGenerators/LayeredVolume.cpp b/MeshLib/MeshGenerators/LayeredVolume.cpp index 569b547ae043286b002cd16e04bf9517da403b8b..254f5b89aff2c0260db3f86063cc23bb42ffe922 100644 --- a/MeshLib/MeshGenerators/LayeredVolume.cpp +++ b/MeshLib/MeshGenerators/LayeredVolume.cpp @@ -15,20 +15,19 @@ #include "LayeredVolume.h" #include "GeoLib/Raster.h" - -#include "MeshLib/Elements/Tri.h" #include "MeshLib/Elements/Quad.h" -#include "MeshLib/Properties.h" +#include "MeshLib/Elements/Tri.h" #include "MeshLib/MeshEditing/DuplicateMeshComponents.h" #include "MeshLib/MeshEditing/RemoveMeshComponents.h" #include "MeshLib/MeshGenerators/MeshLayerMapper.h" #include "MeshLib/MeshSearch/ElementSearch.h" +#include "MeshLib/Properties.h" - -bool LayeredVolume::createRasterLayers(const MeshLib::Mesh &mesh, - const std::vector<GeoLib::Raster const*> &rasters, - double minimum_thickness, - double noDataReplacementValue) +bool LayeredVolume::createRasterLayers( + const MeshLib::Mesh& mesh, + const std::vector<GeoLib::Raster const*>& rasters, + double minimum_thickness, + double noDataReplacementValue) { if (mesh.getDimension() != 2) { @@ -74,7 +73,7 @@ bool LayeredVolume::createRasterLayers(const MeshLib::Mesh &mesh, _materials.resize(_elements.size(), 0); // map each layer and attach to subsurface mesh - const std::size_t nRasters (rasters.size()); + const std::size_t nRasters(rasters.size()); for (std::size_t i = 1; i < nRasters; ++i) { this->addLayerToMesh(*top, i, *rasters[i]); @@ -87,12 +86,14 @@ bool LayeredVolume::createRasterLayers(const MeshLib::Mesh &mesh, return true; } -void LayeredVolume::addLayerToMesh(const MeshLib::Mesh &dem_mesh, unsigned layer_id, GeoLib::Raster const& raster) +void LayeredVolume::addLayerToMesh(const MeshLib::Mesh& dem_mesh, + unsigned layer_id, + GeoLib::Raster const& raster) { - const std::size_t nNodes (dem_mesh.getNumberOfNodes()); - const std::vector<MeshLib::Node*> &nodes (dem_mesh.getNodes()); - const std::size_t node_id_offset (_nodes.size()); - const std::size_t last_layer_node_offset (node_id_offset-nNodes); + const std::size_t nNodes(dem_mesh.getNumberOfNodes()); + const std::vector<MeshLib::Node*>& nodes(dem_mesh.getNodes()); + const std::size_t node_id_offset(_nodes.size()); + const std::size_t last_layer_node_offset(node_id_offset - nNodes); for (std::size_t i = 0; i < nNodes; ++i) { @@ -102,48 +103,57 @@ void LayeredVolume::addLayerToMesh(const MeshLib::Mesh &dem_mesh, unsigned layer _nodes.size())); } - const std::vector<MeshLib::Element*> &layer_elements (dem_mesh.getElements()); + const std::vector<MeshLib::Element*>& layer_elements( + dem_mesh.getElements()); for (MeshLib::Element* elem : layer_elements) { if (elem->getGeomType() == MeshLib::MeshElemType::TRIANGLE) { - std::array<MeshLib::Node*,3> tri_nodes = {{ _nodes[node_id_offset+elem->getNodeIndex(0)], - _nodes[node_id_offset+elem->getNodeIndex(1)], - _nodes[node_id_offset+elem->getNodeIndex(2)] }}; + std::array<MeshLib::Node*, 3> tri_nodes = { + {_nodes[node_id_offset + elem->getNodeIndex(0)], + _nodes[node_id_offset + elem->getNodeIndex(1)], + _nodes[node_id_offset + elem->getNodeIndex(2)]}}; _elements.push_back(new MeshLib::Tri(tri_nodes)); _materials.push_back(layer_id); } else if (elem->getGeomType() == MeshLib::MeshElemType::QUAD) { - std::array<MeshLib::Node*,4> quad_nodes = {{ _nodes[node_id_offset+elem->getNodeIndex(0)], - _nodes[node_id_offset+elem->getNodeIndex(1)], - _nodes[node_id_offset+elem->getNodeIndex(2)], - _nodes[node_id_offset+elem->getNodeIndex(3)] }}; + std::array<MeshLib::Node*, 4> quad_nodes = { + {_nodes[node_id_offset + elem->getNodeIndex(0)], + _nodes[node_id_offset + elem->getNodeIndex(1)], + _nodes[node_id_offset + elem->getNodeIndex(2)], + _nodes[node_id_offset + elem->getNodeIndex(3)]}}; _elements.push_back(new MeshLib::Quad(quad_nodes)); _materials.push_back(layer_id); } } } -void LayeredVolume::addLayerBoundaries(const MeshLib::Mesh &layer, std::size_t nLayers) +void LayeredVolume::addLayerBoundaries(const MeshLib::Mesh& layer, + std::size_t nLayers) { - const unsigned nLayerBoundaries (nLayers-1); - const std::size_t nNodes (layer.getNumberOfNodes()); - const std::vector<MeshLib::Element*> &layer_elements (layer.getElements()); + const unsigned nLayerBoundaries(nLayers - 1); + const std::size_t nNodes(layer.getNumberOfNodes()); + const std::vector<MeshLib::Element*>& layer_elements(layer.getElements()); for (MeshLib::Element* elem : layer_elements) { - const std::size_t nElemNodes (elem->getNumberOfBaseNodes()); + const std::size_t nElemNodes(elem->getNumberOfBaseNodes()); for (unsigned i = 0; i < nElemNodes; ++i) { if (elem->getNeighbor(i) == nullptr) { - for (unsigned j=0; j<nLayerBoundaries; ++j) + for (unsigned j = 0; j < nLayerBoundaries; ++j) { - const std::size_t offset (j*nNodes); + const std::size_t offset(j * nNodes); MeshLib::Node* n0 = _nodes[offset + elem->getNodeIndex(i)]; - MeshLib::Node* n1 = _nodes[offset + elem->getNodeIndex((i+1)%nElemNodes)]; - MeshLib::Node* n2 = _nodes[offset + nNodes + elem->getNodeIndex((i+1)%nElemNodes)]; - MeshLib::Node* n3 = _nodes[offset + nNodes + elem->getNodeIndex(i)]; + MeshLib::Node* n1 = + _nodes[offset + + elem->getNodeIndex((i + 1) % nElemNodes)]; + MeshLib::Node* n2 = + _nodes[offset + nNodes + + elem->getNodeIndex((i + 1) % nElemNodes)]; + MeshLib::Node* n3 = + _nodes[offset + nNodes + elem->getNodeIndex(i)]; auto const v0 = Eigen::Map<Eigen::Vector3d const>(n0->getCoords()); @@ -159,13 +169,13 @@ void LayeredVolume::addLayerBoundaries(const MeshLib::Mesh &layer, std::size_t n { const std::array tri_nodes = {n0, n2, n1}; _elements.push_back(new MeshLib::Tri(tri_nodes)); - _materials.push_back(nLayers+j); + _materials.push_back(nLayers + j); } if ((v3 - v0).norm() > eps) { const std::array tri_nodes = {n0, n3, n2}; _elements.push_back(new MeshLib::Tri(tri_nodes)); - _materials.push_back(nLayers+j); + _materials.push_back(nLayers + j); } } } @@ -173,19 +183,20 @@ void LayeredVolume::addLayerBoundaries(const MeshLib::Mesh &layer, std::size_t n } } -void LayeredVolume::removeCongruentElements(std::size_t nLayers, std::size_t nElementsPerLayer) +void LayeredVolume::removeCongruentElements(std::size_t nLayers, + std::size_t nElementsPerLayer) { - for (std::size_t i=nLayers-1; i>0; --i) + for (std::size_t i = nLayers - 1; i > 0; --i) { - const std::size_t lower_offset ((i-1) * nElementsPerLayer); - const std::size_t upper_offset ( i * nElementsPerLayer); - for (std::size_t j=0; j<nElementsPerLayer; ++j) + const std::size_t lower_offset((i - 1) * nElementsPerLayer); + const std::size_t upper_offset(i * nElementsPerLayer); + for (std::size_t j = 0; j < nElementsPerLayer; ++j) { - MeshLib::Element const*const high (_elements[upper_offset+j]); - MeshLib::Element *const low (_elements[lower_offset+j]); + MeshLib::Element const* const high(_elements[upper_offset + j]); + MeshLib::Element* const low(_elements[lower_offset + j]); unsigned count(0); - const std::size_t nElemNodes (low->getNumberOfBaseNodes()); + const std::size_t nElemNodes(low->getNumberOfBaseNodes()); for (std::size_t k = 0; k < nElemNodes; ++k) { if (high->getNodeIndex(k) == low->getNodeIndex(k)) @@ -197,10 +208,10 @@ void LayeredVolume::removeCongruentElements(std::size_t nLayers, std::size_t nEl if (count == nElemNodes) { - delete _elements[upper_offset+j]; + delete _elements[upper_offset + j]; // mark element and material entries for deletion - _elements[upper_offset+j] = nullptr; - _materials[upper_offset+j] = -1; + _elements[upper_offset + j] = nullptr; + _materials[upper_offset + j] = -1; } else { @@ -214,7 +225,8 @@ void LayeredVolume::removeCongruentElements(std::size_t nLayers, std::size_t nEl } } // delete marked entries - auto elem_vec_end = std::remove(_elements.begin(), _elements.end(), nullptr); + auto elem_vec_end = + std::remove(_elements.begin(), _elements.end(), nullptr); _elements.erase(elem_vec_end, _elements.end()); auto mat_vec_end = std::remove(_materials.begin(), _materials.end(), -1); _materials.erase(mat_vec_end, _materials.end()); diff --git a/MeshLib/MeshGenerators/MeshGenerator.cpp b/MeshLib/MeshGenerators/MeshGenerator.cpp index 6143f2c8633d514769197f84cf642ce4a60acdb9..b85b549b977168427d208fd5c67e7df207f82cc1 100644 --- a/MeshLib/MeshGenerators/MeshGenerator.cpp +++ b/MeshLib/MeshGenerators/MeshGenerator.cpp @@ -25,9 +25,8 @@ namespace MeshLib { - std::vector<MeshLib::Node*> MeshGenerator::generateRegularNodes( - const std::vector<const std::vector<double>*> &vec_xyz_coords, + const std::vector<const std::vector<double>*>& vec_xyz_coords, const MathLib::Point3d& origin) { auto const shift_coordinates = [](auto const& in, auto& out, @@ -58,12 +57,11 @@ std::vector<MeshLib::Node*> MeshGenerator::generateRegularNodes( } std::vector<MeshLib::Node*> MeshGenerator::generateRegularNodes( - const std::vector<double> &vec_x_coords, - const MathLib::Point3d& origin) + const std::vector<double>& vec_x_coords, const MathLib::Point3d& origin) { std::vector<const std::vector<double>*> vec_xyz_coords; vec_xyz_coords.push_back(&vec_x_coords); - std::vector<double> dummy(1,0.0); + std::vector<double> dummy(1, 0.0); for (unsigned i = vec_xyz_coords.size() - 1; i < 3u; i++) { vec_xyz_coords.push_back(&dummy); @@ -72,14 +70,14 @@ std::vector<MeshLib::Node*> MeshGenerator::generateRegularNodes( } std::vector<MeshLib::Node*> MeshGenerator::generateRegularNodes( - std::vector<double> &vec_x_coords, - std::vector<double> &vec_y_coords, + std::vector<double>& vec_x_coords, + std::vector<double>& vec_y_coords, const MathLib::Point3d& origin) { std::vector<const std::vector<double>*> vec_xyz_coords; vec_xyz_coords.push_back(&vec_x_coords); vec_xyz_coords.push_back(&vec_y_coords); - std::vector<double> dummy(1,0.0); + std::vector<double> dummy(1, 0.0); for (unsigned i = vec_xyz_coords.size() - 1; i < 3u; i++) { vec_xyz_coords.push_back(&dummy); @@ -88,9 +86,9 @@ std::vector<MeshLib::Node*> MeshGenerator::generateRegularNodes( } std::vector<MeshLib::Node*> MeshGenerator::generateRegularNodes( - std::vector<double> &vec_x_coords, - std::vector<double> &vec_y_coords, - std::vector<double> &vec_z_coords, + std::vector<double>& vec_x_coords, + std::vector<double>& vec_y_coords, + std::vector<double>& vec_z_coords, const MathLib::Point3d& origin) { std::vector<const std::vector<double>*> vec_xyz_coords; @@ -101,22 +99,22 @@ std::vector<MeshLib::Node*> MeshGenerator::generateRegularNodes( } std::vector<MeshLib::Node*> MeshGenerator::generateRegularNodes( - const std::array<unsigned,3> &n_cells, - const std::array<double,3> &cell_size, + const std::array<unsigned, 3>& n_cells, + const std::array<double, 3>& cell_size, const MathLib::Point3d& origin) { std::vector<Node*> nodes; - nodes.reserve((n_cells[0]+1)*(n_cells[1]+1)*(n_cells[2]+1)); + nodes.reserve((n_cells[0] + 1) * (n_cells[1] + 1) * (n_cells[2] + 1)); - for (std::size_t i = 0; i < n_cells[2]+1; i++) + for (std::size_t i = 0; i < n_cells[2] + 1; i++) { - const double z (origin[2] + cell_size[2] * i); - for (std::size_t j = 0; j < n_cells[1]+1; j++) + const double z(origin[2] + cell_size[2] * i); + for (std::size_t j = 0; j < n_cells[1] + 1; j++) { - const double y (origin[1] + cell_size[1] * j); - for (std::size_t k = 0; k < n_cells[0]+1; k++) + const double y(origin[1] + cell_size[1] * j); + for (std::size_t k = 0; k < n_cells[0] + 1; k++) { - nodes.push_back (new Node(origin[0] + cell_size[0] * k, y, z)); + nodes.push_back(new Node(origin[0] + cell_size[0] * k, y, z)); } } } @@ -156,34 +154,34 @@ std::vector<MeshLib::Node*> generateRegularPyramidTopNodes( } } // end namespace MeshGenerator -Mesh* MeshGenerator::generateLineMesh( - const double length, - const std::size_t subdivision, - const MathLib::Point3d& origin, - std::string const& mesh_name) +Mesh* MeshGenerator::generateLineMesh(const double length, + const std::size_t subdivision, + const MathLib::Point3d& origin, + std::string const& mesh_name) { - return generateLineMesh(subdivision, length/subdivision, origin, mesh_name); + return generateLineMesh(subdivision, length / subdivision, origin, + mesh_name); } -Mesh* MeshGenerator::generateLineMesh( - const unsigned n_cells, - const double cell_size, - MathLib::Point3d const& origin, - std::string const& mesh_name) +Mesh* MeshGenerator::generateLineMesh(const unsigned n_cells, + const double cell_size, + MathLib::Point3d const& origin, + std::string const& mesh_name) { - return generateLineMesh(BaseLib::UniformSubdivision(n_cells*cell_size, n_cells), origin, mesh_name); + return generateLineMesh( + BaseLib::UniformSubdivision(n_cells * cell_size, n_cells), origin, + mesh_name); } -Mesh* MeshGenerator::generateLineMesh( - const BaseLib::ISubdivision &div, - MathLib::Point3d const& origin, - std::string const& mesh_name) +Mesh* MeshGenerator::generateLineMesh(const BaseLib::ISubdivision& div, + MathLib::Point3d const& origin, + std::string const& mesh_name) { const std::vector<double> vec_x(div()); std::vector<Node*> nodes(generateRegularNodes(vec_x, origin)); - //elements - const std::size_t n_cells = nodes.size()-1; + // elements + const std::size_t n_cells = nodes.size() - 1; std::vector<Element*> elements; elements.reserve(n_cells); @@ -195,65 +193,65 @@ Mesh* MeshGenerator::generateLineMesh( return new Mesh(mesh_name, nodes, elements); } -Mesh* MeshGenerator::generateRegularQuadMesh( - const double length, - const std::size_t subdivision, - const MathLib::Point3d& origin, - std::string const& mesh_name) +Mesh* MeshGenerator::generateRegularQuadMesh(const double length, + const std::size_t subdivision, + const MathLib::Point3d& origin, + std::string const& mesh_name) { return generateRegularQuadMesh(subdivision, subdivision, - length/subdivision, length/subdivision, origin, mesh_name); + length / subdivision, length / subdivision, + origin, mesh_name); } -Mesh* MeshGenerator::generateRegularQuadMesh( - const double x_length, - const double y_length, - const std::size_t x_subdivision, - const std::size_t y_subdivision, - const MathLib::Point3d& origin, - std::string const& mesh_name) +Mesh* MeshGenerator::generateRegularQuadMesh(const double x_length, + const double y_length, + const std::size_t x_subdivision, + const std::size_t y_subdivision, + const MathLib::Point3d& origin, + std::string const& mesh_name) { return generateRegularQuadMesh(x_subdivision, y_subdivision, - x_length/x_subdivision, y_length/y_subdivision, origin, mesh_name); + x_length / x_subdivision, + y_length / y_subdivision, origin, mesh_name); } -Mesh* MeshGenerator::generateRegularQuadMesh( - const unsigned n_x_cells, - const unsigned n_y_cells, - const double cell_size, - MathLib::Point3d const& origin, - std::string const& mesh_name) +Mesh* MeshGenerator::generateRegularQuadMesh(const unsigned n_x_cells, + const unsigned n_y_cells, + const double cell_size, + MathLib::Point3d const& origin, + std::string const& mesh_name) { - return generateRegularQuadMesh(n_x_cells, n_y_cells, cell_size, cell_size, origin, mesh_name); + return generateRegularQuadMesh(n_x_cells, n_y_cells, cell_size, cell_size, + origin, mesh_name); } -Mesh* MeshGenerator::generateRegularQuadMesh( - const unsigned n_x_cells, - const unsigned n_y_cells, - const double cell_size_x, - const double cell_size_y, - MathLib::Point3d const& origin, - std::string const& mesh_name) +Mesh* MeshGenerator::generateRegularQuadMesh(const unsigned n_x_cells, + const unsigned n_y_cells, + const double cell_size_x, + const double cell_size_y, + MathLib::Point3d const& origin, + std::string const& mesh_name) { - return generateRegularQuadMesh(BaseLib::UniformSubdivision(n_x_cells*cell_size_x, n_x_cells), - BaseLib::UniformSubdivision(n_y_cells*cell_size_y, n_y_cells), origin, mesh_name); + return generateRegularQuadMesh( + BaseLib::UniformSubdivision(n_x_cells * cell_size_x, n_x_cells), + BaseLib::UniformSubdivision(n_y_cells * cell_size_y, n_y_cells), origin, + mesh_name); } -Mesh* MeshGenerator::generateRegularQuadMesh( - const BaseLib::ISubdivision &div_x, - const BaseLib::ISubdivision &div_y, - MathLib::Point3d const& origin, - std::string const& mesh_name) +Mesh* MeshGenerator::generateRegularQuadMesh(const BaseLib::ISubdivision& div_x, + const BaseLib::ISubdivision& div_y, + MathLib::Point3d const& origin, + std::string const& mesh_name) { std::vector<double> vec_x(div_x()); std::vector<double> vec_y(div_y()); std::vector<Node*> nodes(generateRegularNodes(vec_x, vec_y, origin)); - const unsigned n_x_nodes (vec_x.size()); + const unsigned n_x_nodes(vec_x.size()); - //elements + // elements std::vector<Element*> elements; - const unsigned n_x_cells (vec_x.size()-1); - const unsigned n_y_cells (vec_y.size()-1); + const unsigned n_x_cells(vec_x.size() - 1); + const unsigned n_y_cells(vec_y.size() - 1); elements.reserve(n_x_cells * n_y_cells); for (std::size_t j = 0; j < n_y_cells; j++) @@ -271,85 +269,83 @@ Mesh* MeshGenerator::generateRegularQuadMesh( return new Mesh(mesh_name, nodes, elements); } -Mesh* MeshGenerator::generateRegularHexMesh( - const double length, - const std::size_t subdivision, - const MathLib::Point3d& origin, - std::string const& mesh_name) +Mesh* MeshGenerator::generateRegularHexMesh(const double length, + const std::size_t subdivision, + const MathLib::Point3d& origin, + std::string const& mesh_name) { - return MeshGenerator::generateRegularHexMesh(subdivision, subdivision, - subdivision, length/subdivision, origin, mesh_name); + return MeshGenerator::generateRegularHexMesh( + subdivision, subdivision, subdivision, length / subdivision, origin, + mesh_name); } -Mesh* MeshGenerator::generateRegularHexMesh( - const double x_length, - const double y_length, - const double z_length, - const std::size_t x_subdivision, - const std::size_t y_subdivision, - const std::size_t z_subdivision, - const MathLib::Point3d& origin, - std::string const& mesh_name) +Mesh* MeshGenerator::generateRegularHexMesh(const double x_length, + const double y_length, + const double z_length, + const std::size_t x_subdivision, + const std::size_t y_subdivision, + const std::size_t z_subdivision, + const MathLib::Point3d& origin, + std::string const& mesh_name) { - return MeshGenerator::generateRegularHexMesh(x_subdivision, y_subdivision, z_subdivision, - x_length/x_subdivision, y_length/y_subdivision, z_length/z_subdivision, origin, mesh_name); + return MeshGenerator::generateRegularHexMesh( + x_subdivision, y_subdivision, z_subdivision, x_length / x_subdivision, + y_length / y_subdivision, z_length / z_subdivision, origin, mesh_name); } -Mesh* MeshGenerator::generateRegularHexMesh( - const unsigned n_x_cells, - const unsigned n_y_cells, - const unsigned n_z_cells, - const double cell_size, - MathLib::Point3d const& origin, - std::string const& mesh_name) +Mesh* MeshGenerator::generateRegularHexMesh(const unsigned n_x_cells, + const unsigned n_y_cells, + const unsigned n_z_cells, + const double cell_size, + MathLib::Point3d const& origin, + std::string const& mesh_name) { - return MeshGenerator::generateRegularHexMesh(n_x_cells, n_y_cells, n_z_cells, - cell_size, cell_size, cell_size, origin, mesh_name); + return MeshGenerator::generateRegularHexMesh( + n_x_cells, n_y_cells, n_z_cells, cell_size, cell_size, cell_size, + origin, mesh_name); } -Mesh* MeshGenerator::generateRegularHexMesh( - const unsigned n_x_cells, - const unsigned n_y_cells, - const unsigned n_z_cells, - const double cell_size_x, - const double cell_size_y, - const double cell_size_z, - MathLib::Point3d const& origin, - std::string const& mesh_name) +Mesh* MeshGenerator::generateRegularHexMesh(const unsigned n_x_cells, + const unsigned n_y_cells, + const unsigned n_z_cells, + const double cell_size_x, + const double cell_size_y, + const double cell_size_z, + MathLib::Point3d const& origin, + std::string const& mesh_name) { return generateRegularHexMesh( - BaseLib::UniformSubdivision(n_x_cells*cell_size_x, n_x_cells), - BaseLib::UniformSubdivision(n_y_cells*cell_size_y, n_y_cells), - BaseLib::UniformSubdivision(n_z_cells*cell_size_z, n_z_cells), - origin, mesh_name); + BaseLib::UniformSubdivision(n_x_cells * cell_size_x, n_x_cells), + BaseLib::UniformSubdivision(n_y_cells * cell_size_y, n_y_cells), + BaseLib::UniformSubdivision(n_z_cells * cell_size_z, n_z_cells), origin, + mesh_name); } -Mesh* MeshGenerator::generateRegularHexMesh( - const BaseLib::ISubdivision &div_x, - const BaseLib::ISubdivision &div_y, - const BaseLib::ISubdivision &div_z, - MathLib::Point3d const& origin, - std::string const& mesh_name) +Mesh* MeshGenerator::generateRegularHexMesh(const BaseLib::ISubdivision& div_x, + const BaseLib::ISubdivision& div_y, + const BaseLib::ISubdivision& div_z, + MathLib::Point3d const& origin, + std::string const& mesh_name) { std::vector<double> vec_x(div_x()); std::vector<double> vec_y(div_y()); std::vector<double> vec_z(div_z()); std::vector<Node*> nodes(generateRegularNodes(vec_x, vec_y, vec_z, origin)); - const unsigned n_x_nodes (vec_x.size()); - const unsigned n_y_nodes (vec_y.size()); - const unsigned n_x_cells (vec_x.size()-1); - const unsigned n_y_cells (vec_y.size()-1); - const unsigned n_z_cells (vec_z.size()-1); + const unsigned n_x_nodes(vec_x.size()); + const unsigned n_y_nodes(vec_y.size()); + const unsigned n_x_cells(vec_x.size() - 1); + const unsigned n_y_cells(vec_y.size() - 1); + const unsigned n_z_cells(vec_z.size() - 1); - //elements + // elements std::vector<Element*> elements; elements.reserve(n_x_cells * n_y_cells * n_z_cells); for (std::size_t i = 0; i < n_z_cells; i++) { - const std::size_t offset_z1 = i * n_x_nodes * n_y_nodes; // bottom - const std::size_t offset_z2 = (i + 1) * n_x_nodes * n_y_nodes; // top + const std::size_t offset_z1 = i * n_x_nodes * n_y_nodes; // bottom + const std::size_t offset_z2 = (i + 1) * n_x_nodes * n_y_nodes; // top for (std::size_t j = 0; j < n_y_cells; j++) { const std::size_t offset_y1 = j * n_x_nodes; @@ -470,62 +466,63 @@ Mesh* MeshGenerator::generateRegularPyramidMesh( return new Mesh(mesh_name, nodes, elements); } -Mesh* MeshGenerator::generateRegularTriMesh( - const double length, - const std::size_t subdivision, - const MathLib::Point3d& origin, - std::string const& mesh_name) +Mesh* MeshGenerator::generateRegularTriMesh(const double length, + const std::size_t subdivision, + const MathLib::Point3d& origin, + std::string const& mesh_name) { - return generateRegularTriMesh(subdivision, subdivision, length/subdivision, origin, mesh_name); + return generateRegularTriMesh(subdivision, subdivision, + length / subdivision, origin, mesh_name); } -Mesh* MeshGenerator::generateRegularTriMesh( - const double x_length, - const double y_length, - const std::size_t x_subdivision, - const std::size_t y_subdivision, - const MathLib::Point3d& origin, - std::string const& mesh_name) +Mesh* MeshGenerator::generateRegularTriMesh(const double x_length, + const double y_length, + const std::size_t x_subdivision, + const std::size_t y_subdivision, + const MathLib::Point3d& origin, + std::string const& mesh_name) { - return generateRegularTriMesh(x_subdivision, y_subdivision, x_length/x_subdivision, y_length/y_subdivision, origin, mesh_name); + return generateRegularTriMesh(x_subdivision, y_subdivision, + x_length / x_subdivision, + y_length / y_subdivision, origin, mesh_name); } -Mesh* MeshGenerator::generateRegularTriMesh( - const unsigned n_x_cells, - const unsigned n_y_cells, - const double cell_size, - MathLib::Point3d const& origin, - std::string const& mesh_name) +Mesh* MeshGenerator::generateRegularTriMesh(const unsigned n_x_cells, + const unsigned n_y_cells, + const double cell_size, + MathLib::Point3d const& origin, + std::string const& mesh_name) { - return generateRegularTriMesh(n_x_cells, n_y_cells, cell_size, cell_size, origin, mesh_name); + return generateRegularTriMesh(n_x_cells, n_y_cells, cell_size, cell_size, + origin, mesh_name); } -Mesh* MeshGenerator::generateRegularTriMesh( - const unsigned n_x_cells, - const unsigned n_y_cells, - const double cell_size_x, - const double cell_size_y, - MathLib::Point3d const& origin, - std::string const& mesh_name) +Mesh* MeshGenerator::generateRegularTriMesh(const unsigned n_x_cells, + const unsigned n_y_cells, + const double cell_size_x, + const double cell_size_y, + MathLib::Point3d const& origin, + std::string const& mesh_name) { - return generateRegularTriMesh(BaseLib::UniformSubdivision(n_x_cells*cell_size_x, n_x_cells), - BaseLib::UniformSubdivision(n_y_cells*cell_size_y, n_y_cells), origin, mesh_name); + return generateRegularTriMesh( + BaseLib::UniformSubdivision(n_x_cells * cell_size_x, n_x_cells), + BaseLib::UniformSubdivision(n_y_cells * cell_size_y, n_y_cells), origin, + mesh_name); } -Mesh* MeshGenerator::generateRegularTriMesh( - const BaseLib::ISubdivision &div_x, - const BaseLib::ISubdivision &div_y, - MathLib::Point3d const& origin, - std::string const& mesh_name) +Mesh* MeshGenerator::generateRegularTriMesh(const BaseLib::ISubdivision& div_x, + const BaseLib::ISubdivision& div_y, + MathLib::Point3d const& origin, + std::string const& mesh_name) { std::vector<double> vec_x(div_x()); std::vector<double> vec_y(div_y()); std::vector<Node*> nodes(generateRegularNodes(vec_x, vec_y, origin)); - const unsigned n_x_nodes (vec_x.size()); - const unsigned n_x_cells (vec_x.size()-1); - const unsigned n_y_cells (vec_y.size()-1); + const unsigned n_x_nodes(vec_x.size()); + const unsigned n_x_cells(vec_x.size() - 1); + const unsigned n_y_cells(vec_y.size() - 1); - //elements + // elements std::vector<Element*> elements; elements.reserve(n_x_cells * n_y_cells * 2); @@ -548,115 +545,109 @@ Mesh* MeshGenerator::generateRegularTriMesh( return new Mesh(mesh_name, nodes, elements); } -Mesh* MeshGenerator::generateRegularPrismMesh( - const double x_length, - const double y_length, - const double z_length, - const std::size_t x_subdivision, - const std::size_t y_subdivision, - const std::size_t z_subdivision, - const MathLib::Point3d& origin, - std::string const& mesh_name) +Mesh* MeshGenerator::generateRegularPrismMesh(const double x_length, + const double y_length, + const double z_length, + const std::size_t x_subdivision, + const std::size_t y_subdivision, + const std::size_t z_subdivision, + const MathLib::Point3d& origin, + std::string const& mesh_name) { - return generateRegularPrismMesh(x_subdivision, y_subdivision, z_subdivision, - x_length/x_subdivision, y_length/y_subdivision, z_length/z_subdivision, - origin, mesh_name); + return generateRegularPrismMesh( + x_subdivision, y_subdivision, z_subdivision, x_length / x_subdivision, + y_length / y_subdivision, z_length / z_subdivision, origin, mesh_name); } -Mesh* MeshGenerator::generateRegularPrismMesh( - const unsigned n_x_cells, - const unsigned n_y_cells, - const unsigned n_z_cells, - const double cell_size, - MathLib::Point3d const& origin, - std::string const& mesh_name) +Mesh* MeshGenerator::generateRegularPrismMesh(const unsigned n_x_cells, + const unsigned n_y_cells, + const unsigned n_z_cells, + const double cell_size, + MathLib::Point3d const& origin, + std::string const& mesh_name) { - return generateRegularPrismMesh(n_x_cells, n_y_cells, n_z_cells, - cell_size, cell_size, cell_size, origin, mesh_name); + return generateRegularPrismMesh(n_x_cells, n_y_cells, n_z_cells, cell_size, + cell_size, cell_size, origin, mesh_name); } -Mesh* MeshGenerator::generateRegularPrismMesh( - const unsigned n_x_cells, - const unsigned n_y_cells, - const unsigned n_z_cells, - const double cell_size_x, - const double cell_size_y, - const double cell_size_z, - MathLib::Point3d const& origin, - std::string const& mesh_name) +Mesh* MeshGenerator::generateRegularPrismMesh(const unsigned n_x_cells, + const unsigned n_y_cells, + const unsigned n_z_cells, + const double cell_size_x, + const double cell_size_y, + const double cell_size_z, + MathLib::Point3d const& origin, + std::string const& mesh_name) { - std::unique_ptr<MeshLib::Mesh> mesh ( - generateRegularTriMesh(n_x_cells, n_y_cells, cell_size_x, cell_size_y, origin, mesh_name)); - std::size_t const n_tris (mesh->getNumberOfElements()); + std::unique_ptr<MeshLib::Mesh> mesh(generateRegularTriMesh( + n_x_cells, n_y_cells, cell_size_x, cell_size_y, origin, mesh_name)); + std::size_t const n_tris(mesh->getNumberOfElements()); bool const copy_material_ids = false; for (std::size_t i = 0; i < n_z_cells; ++i) { mesh.reset(MeshLib::addLayerToMesh(*mesh, cell_size_z, mesh_name, true, copy_material_ids)); } - std::vector<std::size_t> elem_ids (n_tris); + std::vector<std::size_t> elem_ids(n_tris); std::iota(elem_ids.begin(), elem_ids.end(), 0); return MeshLib::removeElements(*mesh, elem_ids, mesh_name); } -Mesh* MeshGenerator::generateRegularTetMesh( - const double x_length, - const double y_length, - const double z_length, - const std::size_t x_subdivision, - const std::size_t y_subdivision, - const std::size_t z_subdivision, - const MathLib::Point3d& origin, - std::string const& mesh_name) +Mesh* MeshGenerator::generateRegularTetMesh(const double x_length, + const double y_length, + const double z_length, + const std::size_t x_subdivision, + const std::size_t y_subdivision, + const std::size_t z_subdivision, + const MathLib::Point3d& origin, + std::string const& mesh_name) { - return generateRegularTetMesh(x_subdivision, y_subdivision, z_subdivision, - x_length/x_subdivision, y_length/y_subdivision, z_length/z_subdivision, - origin, mesh_name); + return generateRegularTetMesh( + x_subdivision, y_subdivision, z_subdivision, x_length / x_subdivision, + y_length / y_subdivision, z_length / z_subdivision, origin, mesh_name); } -Mesh* MeshGenerator::generateRegularTetMesh( - const unsigned n_x_cells, - const unsigned n_y_cells, - const unsigned n_z_cells, - const double cell_size_x, - const double cell_size_y, - const double cell_size_z, - MathLib::Point3d const& origin, - std::string const& mesh_name) +Mesh* MeshGenerator::generateRegularTetMesh(const unsigned n_x_cells, + const unsigned n_y_cells, + const unsigned n_z_cells, + const double cell_size_x, + const double cell_size_y, + const double cell_size_z, + MathLib::Point3d const& origin, + std::string const& mesh_name) { return generateRegularTetMesh( - BaseLib::UniformSubdivision(n_x_cells*cell_size_x, n_x_cells), - BaseLib::UniformSubdivision(n_y_cells*cell_size_y, n_y_cells), - BaseLib::UniformSubdivision(n_z_cells*cell_size_z, n_z_cells), - origin, mesh_name); + BaseLib::UniformSubdivision(n_x_cells * cell_size_x, n_x_cells), + BaseLib::UniformSubdivision(n_y_cells * cell_size_y, n_y_cells), + BaseLib::UniformSubdivision(n_z_cells * cell_size_z, n_z_cells), origin, + mesh_name); } -Mesh* MeshGenerator::generateRegularTetMesh( - const BaseLib::ISubdivision &div_x, - const BaseLib::ISubdivision &div_y, - const BaseLib::ISubdivision &div_z, - MathLib::Point3d const& origin, - std::string const& mesh_name) +Mesh* MeshGenerator::generateRegularTetMesh(const BaseLib::ISubdivision& div_x, + const BaseLib::ISubdivision& div_y, + const BaseLib::ISubdivision& div_z, + MathLib::Point3d const& origin, + std::string const& mesh_name) { std::vector<double> vec_x(div_x()); std::vector<double> vec_y(div_y()); std::vector<double> vec_z(div_z()); std::vector<Node*> nodes(generateRegularNodes(vec_x, vec_y, vec_z, origin)); - const unsigned n_x_nodes (vec_x.size()); - const unsigned n_y_nodes (vec_y.size()); - const unsigned n_x_cells (vec_x.size()-1); - const unsigned n_y_cells (vec_y.size()-1); - const unsigned n_z_cells (vec_z.size()-1); + const unsigned n_x_nodes(vec_x.size()); + const unsigned n_y_nodes(vec_y.size()); + const unsigned n_x_cells(vec_x.size() - 1); + const unsigned n_y_cells(vec_y.size() - 1); + const unsigned n_z_cells(vec_z.size() - 1); - //elements + // elements std::vector<Element*> elements; elements.reserve(n_x_cells * n_y_cells * n_z_cells * 6); for (std::size_t i = 0; i < n_z_cells; i++) { - const std::size_t offset_z1 = i * n_x_nodes * n_y_nodes; // bottom - const std::size_t offset_z2 = (i + 1) * n_x_nodes * n_y_nodes; // top + const std::size_t offset_z1 = i * n_x_nodes * n_y_nodes; // bottom + const std::size_t offset_z2 = (i + 1) * n_x_nodes * n_y_nodes; // top for (std::size_t j = 0; j < n_y_cells; j++) { const std::size_t offset_y1 = j * n_x_nodes; @@ -718,18 +709,19 @@ Mesh* MeshGenerator::generateRegularTetMesh( return new Mesh(mesh_name, nodes, elements); } -MeshLib::Mesh* -MeshGenerator::createSurfaceMesh(std::string const& mesh_name, - MathLib::Point3d const& ll, MathLib::Point3d const& ur, - std::array<std::size_t, 2> const& n_steps, - const std::function<double(double,double)>& f) +MeshLib::Mesh* MeshGenerator::createSurfaceMesh( + std::string const& mesh_name, MathLib::Point3d const& ll, + MathLib::Point3d const& ur, std::array<std::size_t, 2> const& n_steps, + const std::function<double(double, double)>& f) { - std::array<double, 2> step_size{{ - (ur[0]-ll[0])/(n_steps[0]-1), (ur[1]-ll[1])/(n_steps[1]-1)}}; + std::array<double, 2> step_size{{(ur[0] - ll[0]) / (n_steps[0] - 1), + (ur[1] - ll[1]) / (n_steps[1] - 1)}}; std::vector<MeshLib::Node*> nodes; - for (std::size_t j(0); j<n_steps[1]; ++j) { - for (std::size_t i(0); i<n_steps[0]; ++i) { + for (std::size_t j(0); j < n_steps[1]; ++j) + { + for (std::size_t i(0); i < n_steps[0]; ++i) + { std::size_t const id = i + j * n_steps[1]; double const x = ll[0] + i * step_size[0]; double const y = ll[1] + j * step_size[1]; @@ -739,12 +731,14 @@ MeshGenerator::createSurfaceMesh(std::string const& mesh_name, } std::vector<MeshLib::Element*> sfc_eles; - for (std::size_t j(0); j<n_steps[1]-1; ++j) { - for (std::size_t i(0); i<n_steps[0]-1; ++i) { - std::size_t id_ll(i+j*n_steps[0]); - std::size_t id_lr(i+1+j*n_steps[0]); - std::size_t id_ul(i+(j+1)*n_steps[0]); - std::size_t id_ur(i+1+(j+1)*n_steps[0]); + for (std::size_t j(0); j < n_steps[1] - 1; ++j) + { + for (std::size_t i(0); i < n_steps[0] - 1; ++i) + { + std::size_t id_ll(i + j * n_steps[0]); + std::size_t id_lr(i + 1 + j * n_steps[0]); + std::size_t id_ul(i + (j + 1) * n_steps[0]); + std::size_t id_ur(i + 1 + (j + 1) * n_steps[0]); sfc_eles.push_back( new MeshLib::Tri({nodes[id_ll], nodes[id_lr], nodes[id_ur]})); sfc_eles.push_back( diff --git a/MeshLib/MeshGenerators/MeshLayerMapper.cpp b/MeshLib/MeshGenerators/MeshLayerMapper.cpp index 893af8178d9325a38d7cbf7c7fd4f66615d7e2ec..1c7f87437852e44568239164252a20efcd0298e4 100644 --- a/MeshLib/MeshGenerators/MeshLayerMapper.cpp +++ b/MeshLib/MeshGenerators/MeshLayerMapper.cpp @@ -17,22 +17,20 @@ #include <algorithm> #include "BaseLib/Logging.h" - #include "GeoLib/Raster.h" - #include "MathLib/MathTools.h" - -#include "MeshLib/Elements/Tet.h" #include "MeshLib/Elements/Hex.h" -#include "MeshLib/Elements/Pyramid.h" #include "MeshLib/Elements/Prism.h" +#include "MeshLib/Elements/Pyramid.h" +#include "MeshLib/Elements/Tet.h" #include "MeshLib/MeshSurfaceExtraction.h" #include "MeshLib/Properties.h" namespace MeshLib { - -MeshLib::Mesh* MeshLayerMapper::createStaticLayers(MeshLib::Mesh const& mesh, std::vector<float> const& layer_thickness_vector, std::string const& mesh_name) +MeshLib::Mesh* MeshLayerMapper::createStaticLayers( + MeshLib::Mesh const& mesh, std::vector<float> const& layer_thickness_vector, + std::string const& mesh_name) { std::vector<float> thickness; for (std::size_t i = 0; i < layer_thickness_vector.size(); ++i) @@ -51,18 +49,22 @@ MeshLib::Mesh* MeshLayerMapper::createStaticLayers(MeshLib::Mesh const& mesh, st const std::size_t nLayers(thickness.size()); if (nLayers < 1 || mesh.getDimension() != 2) { - ERR("MeshLayerMapper::createStaticLayers(): A 2D mesh with nLayers > 0 is required as input."); + ERR("MeshLayerMapper::createStaticLayers(): A 2D mesh with nLayers > 0 " + "is required as input."); return nullptr; } const std::size_t nNodes = mesh.getNumberOfNodes(); // count number of 2d elements in the original mesh - const std::size_t nElems (std::count_if(mesh.getElements().begin(), mesh.getElements().end(), - [](MeshLib::Element const* elem) { return (elem->getDimension() == 2);})); - - const std::size_t nOrgElems (mesh.getNumberOfElements()); - const std::vector<MeshLib::Node*> &nodes = mesh.getNodes(); - const std::vector<MeshLib::Element*> &elems = mesh.getElements(); + const std::size_t nElems( + std::count_if(mesh.getElements().begin(), mesh.getElements().end(), + [](MeshLib::Element const* elem) { + return (elem->getDimension() == 2); + })); + + const std::size_t nOrgElems(mesh.getNumberOfElements()); + 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; new_elems.reserve(nElems * nLayers); @@ -76,32 +78,37 @@ MeshLib::Mesh* MeshLayerMapper::createStaticLayers(MeshLib::Mesh const& mesh, st } materials->reserve(nElems * nLayers); - double z_offset (0.0); + 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 node_offset(nNodes * layer_id); if (layer_id > 0) { z_offset += thickness[layer_id - 1]; } - std::transform(nodes.cbegin(), nodes.cend(), new_nodes.begin() + node_offset, - [&z_offset](MeshLib::Node* node){ return new MeshLib::Node((*node)[0], (*node)[1], (*node)[2]-z_offset); }); + std::transform(nodes.cbegin(), nodes.cend(), + new_nodes.begin() + node_offset, + [&z_offset](MeshLib::Node* node) { + return new MeshLib::Node((*node)[0], (*node)[1], + (*node)[2] - z_offset); + }); - // starting with 2nd layer create prism or hex elements connecting the last layer with the current one + // starting with 2nd layer create prism or hex elements connecting the + // last layer with the current one if (layer_id == 0) { continue; } node_offset -= nNodes; - const unsigned mat_id (nLayers - layer_id); + const unsigned mat_id(nLayers - layer_id); for (unsigned i = 0; i < nOrgElems; ++i) { - const MeshLib::Element* sfc_elem( elems[i] ); + const MeshLib::Element* sfc_elem(elems[i]); if (sfc_elem->getDimension() < 2) { // ignore line-elements continue; @@ -110,11 +117,12 @@ MeshLib::Mesh* MeshLayerMapper::createStaticLayers(MeshLib::Mesh const& mesh, st const unsigned nElemNodes(sfc_elem->getNumberOfBaseNodes()); auto** e_nodes = new MeshLib::Node*[2 * nElemNodes]; - for (unsigned j=0; j<nElemNodes; ++j) + 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]; + 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() == MeshLib::MeshElemType::TRIANGLE) { @@ -145,7 +153,8 @@ bool MeshLayerMapper::createRasterLayers( const std::size_t nLayers(rasters.size()); if (nLayers < 2 || mesh.getDimension() != 2) { - ERR("MeshLayerMapper::createRasterLayers(): A 2D mesh and at least two rasters required as input."); + ERR("MeshLayerMapper::createRasterLayers(): A 2D mesh and at least two " + "rasters required as input."); return false; } @@ -166,11 +175,13 @@ bool MeshLayerMapper::createRasterLayers( _nodes.reserve(nLayers * nNodes); // number of triangles in the original mesh - std::size_t const nElems (std::count_if(mesh.getElements().begin(), mesh.getElements().end(), - [](MeshLib::Element const* elem) - { return (elem->getGeomType() == MeshLib::MeshElemType::TRIANGLE);})); - _elements.reserve(nElems * (nLayers-1)); - _materials.reserve(nElems * (nLayers-1)); + std::size_t const nElems(std::count_if( + mesh.getElements().begin(), mesh.getElements().end(), + [](MeshLib::Element const* elem) { + return (elem->getGeomType() == MeshLib::MeshElemType::TRIANGLE); + })); + _elements.reserve(nElems * (nLayers - 1)); + _materials.reserve(nElems * (nLayers - 1)); // add bottom layer std::vector<MeshLib::Node*> const& nodes = bottom->getNodes(); @@ -186,13 +197,14 @@ bool MeshLayerMapper::createRasterLayers( return true; } -void MeshLayerMapper::addLayerToMesh(const MeshLib::Mesh &dem_mesh, unsigned layer_id, GeoLib::Raster const& raster) +void MeshLayerMapper::addLayerToMesh(const MeshLib::Mesh& dem_mesh, + unsigned layer_id, + GeoLib::Raster const& raster) { - const unsigned pyramid_base[3][4] = - { - {1, 3, 4, 2}, // Point 4 missing - {2, 4, 3, 0}, // Point 5 missing - {0, 3, 4, 1}, // Point 6 missing + const unsigned pyramid_base[3][4] = { + {1, 3, 4, 2}, // Point 4 missing + {2, 4, 3, 0}, // Point 5 missing + {0, 3, 4, 1}, // Point 6 missing }; std::size_t const nNodes = dem_mesh.getNumberOfNodes(); @@ -208,11 +220,11 @@ void MeshLayerMapper::addLayerToMesh(const MeshLib::Mesh &dem_mesh, unsigned lay } std::vector<MeshLib::Element*> const& elems = dem_mesh.getElements(); - std::size_t const nElems (dem_mesh.getNumberOfElements()); + std::size_t const nElems(dem_mesh.getNumberOfElements()); - for (std::size_t i=0; i<nElems; ++i) + for (std::size_t i = 0; i < nElems; ++i) { - MeshLib::Element* elem (elems[i]); + MeshLib::Element* elem(elems[i]); if (elem->getGeomType() != MeshLib::MeshElemType::TRIANGLE) { continue; @@ -220,10 +232,14 @@ void MeshLayerMapper::addLayerToMesh(const MeshLib::Mesh &dem_mesh, unsigned lay unsigned node_counter(3); unsigned missing_idx(0); std::array<MeshLib::Node*, 6> new_elem_nodes{}; - for (unsigned j=0; j<3; ++j) + for (unsigned j = 0; j < 3; ++j) { - new_elem_nodes[j] = _nodes[_nodes[last_layer_node_offset + elem->getNodeIndex(j)]->getID()]; - new_elem_nodes[node_counter] = (_nodes[last_layer_node_offset + elem->getNodeIndex(j) + nNodes]); + new_elem_nodes[j] = + _nodes[_nodes[last_layer_node_offset + elem->getNodeIndex(j)] + ->getID()]; + new_elem_nodes[node_counter] = + (_nodes[last_layer_node_offset + elem->getNodeIndex(j) + + nNodes]); if (new_elem_nodes[j]->getID() != new_elem_nodes[node_counter]->getID()) { @@ -237,32 +253,34 @@ void MeshLayerMapper::addLayerToMesh(const MeshLib::Mesh &dem_mesh, unsigned lay switch (node_counter) { - case 6: - _elements.push_back(new MeshLib::Prism(new_elem_nodes)); - _materials.push_back(layer_id); - break; - case 5: - { - std::array<MeshLib::Node*, 5> pyramid_nodes; - pyramid_nodes[0] = new_elem_nodes[pyramid_base[missing_idx][0]]; - pyramid_nodes[1] = new_elem_nodes[pyramid_base[missing_idx][1]]; - pyramid_nodes[2] = new_elem_nodes[pyramid_base[missing_idx][2]]; - pyramid_nodes[3] = new_elem_nodes[pyramid_base[missing_idx][3]]; - pyramid_nodes[4] = new_elem_nodes[missing_idx]; - _elements.push_back(new MeshLib::Pyramid(pyramid_nodes)); - _materials.push_back(layer_id); - break; - } - case 4: - { - std::array<MeshLib::Node*, 4> tet_nodes; - std::copy(new_elem_nodes.begin(), new_elem_nodes.begin() + node_counter, tet_nodes.begin()); - _elements.push_back(new MeshLib::Tet(tet_nodes)); - _materials.push_back(layer_id); - break; - } - default: - continue; + case 6: + _elements.push_back(new MeshLib::Prism(new_elem_nodes)); + _materials.push_back(layer_id); + break; + case 5: + { + std::array<MeshLib::Node*, 5> pyramid_nodes; + pyramid_nodes[0] = new_elem_nodes[pyramid_base[missing_idx][0]]; + pyramid_nodes[1] = new_elem_nodes[pyramid_base[missing_idx][1]]; + pyramid_nodes[2] = new_elem_nodes[pyramid_base[missing_idx][2]]; + pyramid_nodes[3] = new_elem_nodes[pyramid_base[missing_idx][3]]; + pyramid_nodes[4] = new_elem_nodes[missing_idx]; + _elements.push_back(new MeshLib::Pyramid(pyramid_nodes)); + _materials.push_back(layer_id); + break; + } + case 4: + { + std::array<MeshLib::Node*, 4> tet_nodes; + std::copy(new_elem_nodes.begin(), + new_elem_nodes.begin() + node_counter, + tet_nodes.begin()); + _elements.push_back(new MeshLib::Tet(tet_nodes)); + _materials.push_back(layer_id); + break; + } + default: + continue; } } } @@ -272,23 +290,24 @@ bool MeshLayerMapper::layerMapping(MeshLib::Mesh const& mesh, double const nodata_replacement, bool const ignore_nodata) { - if (mesh.getDimension() != 2) { ERR("MshLayerMapper::layerMapping() - requires 2D mesh"); return false; } - GeoLib::RasterHeader const& header (raster.getHeader()); + GeoLib::RasterHeader const& header(raster.getHeader()); const double x0(header.origin[0]); const double y0(header.origin[1]); const double delta(header.cell_size); - const std::pair<double, double> xDim(x0, x0 + header.n_cols * delta); // extension in x-dimension - const std::pair<double, double> yDim(y0, y0 + header.n_rows * delta); // extension in y-dimension + const std::pair<double, double> xDim( + x0, x0 + header.n_cols * delta); // extension in x-dimension + const std::pair<double, double> yDim( + y0, y0 + header.n_rows * delta); // extension in y-dimension - const std::size_t nNodes (mesh.getNumberOfNodes()); - const std::vector<MeshLib::Node*> &nodes = mesh.getNodes(); + const std::size_t nNodes(mesh.getNumberOfNodes()); + const std::vector<MeshLib::Node*>& nodes = mesh.getNodes(); for (unsigned i = 0; i < nNodes; ++i) { if (!ignore_nodata && !raster.isPntOnRaster(*nodes[i])) @@ -299,7 +318,7 @@ bool MeshLayerMapper::layerMapping(MeshLib::Mesh const& mesh, continue; } - double elevation (raster.getValueAtPoint(*nodes[i])); + double elevation(raster.getValueAtPoint(*nodes[i])); constexpr double eps = std::numeric_limits<double>::epsilon(); if (std::abs(elevation - header.no_data) < eps) { @@ -328,7 +347,7 @@ bool MeshLayerMapper::mapToStaticValue(MeshLib::Mesh const& mesh, return false; } - std::vector<MeshLib::Node*> const& nodes (mesh.getNodes()); + std::vector<MeshLib::Node*> const& nodes(mesh.getNodes()); for (MeshLib::Node* node : nodes) { node->updateCoordinates((*node)[0], (*node)[1], value); @@ -336,4 +355,4 @@ bool MeshLayerMapper::mapToStaticValue(MeshLib::Mesh const& mesh, return true; } -} // end namespace MeshLib +} // end namespace MeshLib diff --git a/MeshLib/MeshGenerators/RasterToMesh.cpp b/MeshLib/MeshGenerators/RasterToMesh.cpp index e010688f0ee4efe2a8591be7841b3750823a8786..4599e8b32ab917509dab8eed87d95ea9e990e58a 100644 --- a/MeshLib/MeshGenerators/RasterToMesh.cpp +++ b/MeshLib/MeshGenerators/RasterToMesh.cpp @@ -10,21 +10,20 @@ #include "RasterToMesh.h" -#include "BaseLib/Logging.h" +#include <vtkCell.h> +#include <vtkCellData.h> +#include <vtkDataArray.h> +#include <vtkImageData.h> +#include <vtkPointData.h> +#include <vtkSmartPointer.h> +#include "BaseLib/Logging.h" #include "MeshLib/Elements/Elements.h" #include "MeshLib/Mesh.h" -#include "MeshLib/Node.h" #include "MeshLib/MeshEditing/RemoveMeshComponents.h" #include "MeshLib/MeshGenerators/MeshGenerator.h" #include "MeshLib/MeshSearch/ElementSearch.h" - -#include <vtkImageData.h> -#include <vtkDataArray.h> -#include <vtkCell.h> -#include <vtkCellData.h> -#include <vtkSmartPointer.h> -#include <vtkPointData.h> +#include "MeshLib/Node.h" namespace MeshLib { @@ -34,7 +33,11 @@ std::unique_ptr<MeshLib::Mesh> RasterToMesh::convert( UseIntensityAs intensity_type, std::string const& array_name) { - return convert(raster.begin(), raster.getHeader(), elem_type, intensity_type, array_name); + return convert(raster.begin(), + raster.getHeader(), + elem_type, + intensity_type, + array_name); } std::unique_ptr<MeshLib::Mesh> RasterToMesh::convert( @@ -55,10 +58,12 @@ std::unique_ptr<MeshLib::Mesh> RasterToMesh::convert( } int* dims = img->GetDimensions(); - if (((elem_type == MeshElemType::TRIANGLE) || (elem_type == MeshElemType::QUAD)) && + if (((elem_type == MeshElemType::TRIANGLE) || + (elem_type == MeshElemType::QUAD)) && dims[2] != 1) { - ERR("Triangle or Quad elements cannot be used to construct meshes from 3D rasters."); + ERR("Triangle or Quad elements cannot be used to construct meshes from " + "3D rasters."); return nullptr; } @@ -67,7 +72,8 @@ std::unique_ptr<MeshLib::Mesh> RasterToMesh::convert( int nTuple = pixelData->GetNumberOfComponents(); if (nTuple < 1 || nTuple > 4) { - ERR("VtkMeshConverter::convertImgToMesh(): Unsupported pixel composition!"); + ERR("VtkMeshConverter::convertImgToMesh(): Unsupported pixel " + "composition!"); return nullptr; } @@ -84,14 +90,16 @@ std::unique_ptr<MeshLib::Mesh> RasterToMesh::convert( std::vector<double> pix(header.n_cols * header.n_rows * header.n_depth, 0); for (std::size_t k = 0; k < header.n_depth; k++) { - std::size_t const layer_idx = (k*header.n_rows*header.n_cols); + std::size_t const layer_idx = (k * header.n_rows * header.n_cols); for (std::size_t i = 0; i < header.n_rows; i++) { std::size_t const idx = i * header.n_cols + layer_idx; for (std::size_t j = 0; j < header.n_cols; j++) { double* colour = pixelData->GetTuple(idx + j); - bool const visible = (nTuple == 2 || nTuple == 4) ? (colour[nTuple - 1] != 0) : true; + bool const visible = (nTuple == 2 || nTuple == 4) + ? (colour[nTuple - 1] != 0) + : true; if (!visible) { pix[idx + j] = header.no_data; @@ -110,7 +118,7 @@ std::unique_ptr<MeshLib::Mesh> RasterToMesh::convert( } std::unique_ptr<MeshLib::Mesh> RasterToMesh::convert( - double const*const img, + double const* const img, GeoLib::RasterHeader const& header, MeshElemType elem_type, UseIntensityAs intensity_type, @@ -125,15 +133,18 @@ std::unique_ptr<MeshLib::Mesh> RasterToMesh::convert( return nullptr; } - if (((elem_type == MeshElemType::TRIANGLE) || (elem_type == MeshElemType::QUAD)) && + if (((elem_type == MeshElemType::TRIANGLE) || + (elem_type == MeshElemType::QUAD)) && header.n_depth != 1) { - ERR("Triangle or Quad elements cannot be used to construct meshes from 3D rasters."); + ERR("Triangle or Quad elements cannot be used to construct meshes from " + "3D rasters."); return nullptr; } if (intensity_type == UseIntensityAs::ELEVATION && - ((elem_type == MeshElemType::PRISM) || (elem_type == MeshElemType::HEXAHEDRON))) + ((elem_type == MeshElemType::PRISM) || + (elem_type == MeshElemType::HEXAHEDRON))) { ERR("Elevation mapping can only be performed for 2D meshes."); return nullptr; @@ -186,11 +197,12 @@ std::unique_ptr<MeshLib::Mesh> RasterToMesh::convert( std::vector<MeshLib::Node*> const& nodes(mesh->getNodes()); std::vector<MeshLib::Element*> const& elems(mesh->getElements()); std::size_t const n_nodes(elems[0]->getNumberOfNodes()); - bool const double_idx = (elem_type == MeshElemType::TRIANGLE) || (elem_type == MeshElemType::PRISM); + bool const double_idx = (elem_type == MeshElemType::TRIANGLE) || + (elem_type == MeshElemType::PRISM); std::size_t const m = (double_idx) ? 2 : 1; for (std::size_t k = 0; k < header.n_depth; k++) { - std::size_t const layer_idx = (k*header.n_rows*header.n_cols); + std::size_t const layer_idx = (k * header.n_rows * header.n_cols); for (std::size_t i = 0; i < header.n_cols; i++) { std::size_t const idx(i * header.n_rows + layer_idx); @@ -208,7 +220,8 @@ std::unique_ptr<MeshLib::Mesh> RasterToMesh::convert( } for (std::size_t n = 0; n < n_nodes; ++n) { - (*(nodes[elems[m * (idx + j)]->getNodeIndex(n)]))[2] = val; + (*(nodes[elems[m * (idx + j)]->getNodeIndex(n)]))[2] = + val; if (double_idx) { (*(nodes[elems[m * (idx + j) + 1]->getNodeIndex( @@ -221,14 +234,15 @@ std::unique_ptr<MeshLib::Mesh> RasterToMesh::convert( } else { - MeshLib::Properties &properties = mesh->getProperties(); + MeshLib::Properties& properties = mesh->getProperties(); MeshLib::ElementSearch ex(*mesh); if (array_name == "MaterialIDs") { auto* const prop_vec = properties.createNewPropertyVector<int>( array_name, MeshLib::MeshItemType::Cell, 1); fillPropertyVector<int>(*prop_vec, img, header, elem_type); - ex.searchByPropertyValue<int>(array_name, static_cast<int>(header.no_data)); + ex.searchByPropertyValue<int>(array_name, + static_cast<int>(header.no_data)); } else { @@ -257,4 +271,4 @@ std::unique_ptr<MeshLib::Mesh> RasterToMesh::convert( return new_mesh; } -} // end namespace MeshLib +} // end namespace MeshLib diff --git a/MeshLib/MeshGenerators/VtkMeshConverter.cpp b/MeshLib/MeshGenerators/VtkMeshConverter.cpp index 8204821ae995282b13303bc61b528e8951ec90fd..8c97daef1efe0a30350dc01bd2f501bfcd999650 100644 --- a/MeshLib/MeshGenerators/VtkMeshConverter.cpp +++ b/MeshLib/MeshGenerators/VtkMeshConverter.cpp @@ -322,15 +322,14 @@ void VtkMeshConverter::convertArray(vtkDataArray& array, if (vtkLongLongArray::SafeDownCast(&array)) { - VtkMeshConverter::convertTypedArray<long long>( - array, properties, type); + VtkMeshConverter::convertTypedArray<long long>(array, properties, type); return; } if (vtkUnsignedLongArray::SafeDownCast(&array)) { - VtkMeshConverter::convertTypedArray<unsigned long>( - array, properties, type); + VtkMeshConverter::convertTypedArray<unsigned long>(array, properties, + type); return; } diff --git a/MeshLib/MeshQuality/AngleSkewMetric.cpp b/MeshLib/MeshQuality/AngleSkewMetric.cpp index a52b0b0f1c6ef3a2de17c75e56b650a3ed6c6678..a8f77e37c2b5fa7c8b8edd4b77fc774bd5780468 100644 --- a/MeshLib/MeshQuality/AngleSkewMetric.cpp +++ b/MeshLib/MeshQuality/AngleSkewMetric.cpp @@ -14,18 +14,16 @@ #include "AngleSkewMetric.h" -#include <cmath> #include <boost/math/constants/constants.hpp> - -#include "MeshLib/Node.h" +#include <cmath> #include "MathLib/MathTools.h" +#include "MeshLib/Node.h" using namespace boost::math::double_constants; namespace MeshLib { - namespace { template <unsigned long N> @@ -144,9 +142,9 @@ double checkPrism(Element const& elem) } // end unnamed namespace -AngleSkewMetric::AngleSkewMetric(Mesh const& mesh) : - ElementQualityMetric(mesh) -{} +AngleSkewMetric::AngleSkewMetric(Mesh const& mesh) : ElementQualityMetric(mesh) +{ +} void AngleSkewMetric::calculateQuality() { @@ -178,4 +176,4 @@ void AngleSkewMetric::calculateQuality() } } -} // end namespace MeshLib +} // end namespace MeshLib diff --git a/MeshLib/MeshQuality/EdgeRatioMetric.cpp b/MeshLib/MeshQuality/EdgeRatioMetric.cpp index 4e24c4b9c8647d7bbd1fb3ddc852ca8ae0a3ff3b..7ad4a89ba67612d945f435282adc826e7c0a02fa 100644 --- a/MeshLib/MeshQuality/EdgeRatioMetric.cpp +++ b/MeshLib/MeshQuality/EdgeRatioMetric.cpp @@ -19,8 +19,7 @@ namespace MeshLib { -EdgeRatioMetric::EdgeRatioMetric(Mesh const& mesh) : - ElementQualityMetric(mesh) +EdgeRatioMetric::EdgeRatioMetric(Mesh const& mesh) : ElementQualityMetric(mesh) { } @@ -55,4 +54,4 @@ void EdgeRatioMetric::calculateQuality() } } -} // end namespace MeshLib +} // end namespace MeshLib diff --git a/MeshLib/MeshQuality/ElementSizeMetric.cpp b/MeshLib/MeshQuality/ElementSizeMetric.cpp index 29daa1c3731390f5acef61220b75e69a97f59118..3e6cd80225def3b0d375c7e8e8bc08c37055a164 100644 --- a/MeshLib/MeshQuality/ElementSizeMetric.cpp +++ b/MeshLib/MeshQuality/ElementSizeMetric.cpp @@ -19,8 +19,9 @@ namespace MeshLib { ElementSizeMetric::ElementSizeMetric(Mesh const& mesh) -: ElementQualityMetric(mesh) -{} + : ElementQualityMetric(mesh) +{ +} void ElementSizeMetric::calculateQuality() { @@ -50,7 +51,7 @@ void ElementSizeMetric::calculateQuality() std::size_t ElementSizeMetric::calc1dQuality() { - const std::vector<MeshLib::Element*> &elements(_mesh.getElements()); + const std::vector<MeshLib::Element*>& elements(_mesh.getElements()); const std::size_t nElems(elements.size()); std::size_t error_count(0); @@ -79,13 +80,13 @@ std::size_t ElementSizeMetric::calc1dQuality() std::size_t ElementSizeMetric::calc2dQuality() { - const std::vector<MeshLib::Element*> &elements(_mesh.getElements()); + const std::vector<MeshLib::Element*>& elements(_mesh.getElements()); const std::size_t nElems(elements.size()); std::size_t error_count(0); for (std::size_t k(0); k < nElems; k++) { - Element const& elem (*elements[k]); + Element const& elem(*elements[k]); if (elem.getDimension() == 1) { @@ -120,14 +121,14 @@ std::size_t ElementSizeMetric::calc3dQuality() for (std::size_t k(0); k < nElems; k++) { - Element const& elem (*elements[k]); - if (elem.getDimension()<3) + Element const& elem(*elements[k]); + if (elem.getDimension() < 3) { _element_quality_metric[k] = 0.0; continue; } - double const volume (elem.getContent()); + double const volume(elem.getContent()); if (volume < sqrt(std::abs(std::numeric_limits<double>::epsilon()))) { error_count++; @@ -146,4 +147,4 @@ std::size_t ElementSizeMetric::calc3dQuality() return error_count; } -} // end namespace MeshLib +} // end namespace MeshLib diff --git a/MeshLib/MeshQuality/MeshValidation.cpp b/MeshLib/MeshQuality/MeshValidation.cpp index 46b5e64b0b9c94a3dc05cb9e1dbbb5e3656e76d3..daa85d534c708d94f725e77614c4fef6d6a3187c 100644 --- a/MeshLib/MeshQuality/MeshValidation.cpp +++ b/MeshLib/MeshQuality/MeshValidation.cpp @@ -14,28 +14,28 @@ #include "MeshValidation.h" -#include <numeric> #include <algorithm> +#include <numeric> #include <stack> #include "BaseLib/Logging.h" - -#include "MeshLib/Mesh.h" -#include "MeshLib/Node.h" #include "MeshLib/Elements/Element.h" +#include "MeshLib/Mesh.h" #include "MeshLib/MeshEditing/MeshRevision.h" #include "MeshLib/MeshSearch/NodeSearch.h" #include "MeshLib/MeshSurfaceExtraction.h" +#include "MeshLib/Node.h" -namespace MeshLib { - -MeshValidation::MeshValidation(MeshLib::Mesh &mesh) +namespace MeshLib { - INFO ("Mesh Quality Control:"); - INFO ("Looking for unused nodes..."); +MeshValidation::MeshValidation(MeshLib::Mesh& mesh) +{ + INFO("Mesh Quality Control:"); + INFO("Looking for unused nodes..."); NodeSearch ns(mesh); ns.searchUnused(); - if (!ns.getSearchedNodeIDs().empty()) { + if (!ns.getSearchedNodeIDs().empty()) + { INFO("{:d} unused mesh nodes found.", ns.getSearchedNodeIDs().size()); } MeshRevision rev(mesh); @@ -53,19 +53,20 @@ MeshValidation::MeshValidation(MeshLib::Mesh &mesh) } } -std::vector<ElementErrorCode> MeshValidation::testElementGeometry(const MeshLib::Mesh &mesh, double min_volume) +std::vector<ElementErrorCode> MeshValidation::testElementGeometry( + const MeshLib::Mesh& mesh, double min_volume) { - INFO ("Testing mesh element geometry:"); + INFO("Testing mesh element geometry:"); const auto nErrorCodes( static_cast<std::size_t>(ElementErrorFlag::MaxValue)); unsigned error_count[nErrorCodes]; std::fill_n(error_count, 4, 0); - const std::size_t nElements (mesh.getNumberOfElements()); - const std::vector<MeshLib::Element*> &elements (mesh.getElements()); + const std::size_t nElements(mesh.getNumberOfElements()); + const std::vector<MeshLib::Element*>& elements(mesh.getElements()); std::vector<ElementErrorCode> error_code_vector; error_code_vector.reserve(nElements); - for (std::size_t i=0; i<nElements; ++i) + for (std::size_t i = 0; i < nElements; ++i) { const ElementErrorCode e = elements[i]->validate(); error_code_vector.push_back(e); @@ -75,14 +76,17 @@ std::vector<ElementErrorCode> MeshValidation::testElementGeometry(const MeshLib: } // increment error statistics - const std::bitset< static_cast<std::size_t>(ElementErrorFlag::MaxValue) > flags (static_cast< std::bitset<static_cast<std::size_t>(ElementErrorFlag::MaxValue)> >(e)); + const std::bitset<static_cast<std::size_t>(ElementErrorFlag::MaxValue)> + flags(static_cast<std::bitset<static_cast<std::size_t>( + ElementErrorFlag::MaxValue)>>(e)); for (unsigned j = 0; j < nErrorCodes; ++j) { error_count[j] += flags[j]; } } - // if a larger volume threshold is given, evaluate elements again to add them even if they are formally okay + // if a larger volume threshold is given, evaluate elements again to add + // them even if they are formally okay if (min_volume > std::numeric_limits<double>::epsilon()) { for (std::size_t i = 0; i < nElements; ++i) @@ -99,8 +103,9 @@ std::vector<ElementErrorCode> MeshValidation::testElementGeometry(const MeshLib: std::accumulate(error_count, error_count + nErrorCodes, 0.0))); if (error_sum != 0) { - ElementErrorFlag flags[nErrorCodes] = { ElementErrorFlag::ZeroVolume, ElementErrorFlag::NonCoplanar, - ElementErrorFlag::NonConvex, ElementErrorFlag::NodeOrder }; + ElementErrorFlag flags[nErrorCodes] = { + ElementErrorFlag::ZeroVolume, ElementErrorFlag::NonCoplanar, + ElementErrorFlag::NonConvex, ElementErrorFlag::NodeOrder}; for (std::size_t i = 0; i < nErrorCodes; ++i) { if (error_count[i]) @@ -113,13 +118,14 @@ std::vector<ElementErrorCode> MeshValidation::testElementGeometry(const MeshLib: } else { - INFO ("No errors found."); + INFO("No errors found."); } return error_code_vector; } std::array<std::string, static_cast<std::size_t>(ElementErrorFlag::MaxValue)> -MeshValidation::ElementErrorCodeOutput(const std::vector<ElementErrorCode> &error_codes) +MeshValidation::ElementErrorCodeOutput( + const std::vector<ElementErrorCode>& error_codes) { const auto nErrorFlags( static_cast<std::size_t>(ElementErrorFlag::MaxValue)); @@ -188,19 +194,21 @@ unsigned MeshValidation::detectHoles(MeshLib::Mesh const& mesh) return (--current_surface_id); } -void MeshValidation::trackSurface(MeshLib::Element const* element, std::vector<unsigned> &sfc_idx, unsigned const current_index) +void MeshValidation::trackSurface(MeshLib::Element const* element, + std::vector<unsigned>& sfc_idx, + unsigned const current_index) { std::stack<MeshLib::Element const*> elem_stack; elem_stack.push(element); while (!elem_stack.empty()) { - MeshLib::Element const*const elem = elem_stack.top(); + MeshLib::Element const* const elem = elem_stack.top(); elem_stack.pop(); sfc_idx[elem->getID()] = current_index; - std::size_t const n_neighbors (elem->getNumberOfNeighbors()); - for (std::size_t i=0; i<n_neighbors; ++i) + std::size_t const n_neighbors(elem->getNumberOfNeighbors()); + for (std::size_t i = 0; i < n_neighbors; ++i) { - MeshLib::Element const* neighbor (elem->getNeighbor(i)); + MeshLib::Element const* neighbor(elem->getNeighbor(i)); if (neighbor != nullptr && sfc_idx[neighbor->getID()] == std::numeric_limits<unsigned>::max()) { @@ -210,4 +218,4 @@ void MeshValidation::trackSurface(MeshLib::Element const* element, std::vector<u } } -} // end namespace MeshLib +} // end namespace MeshLib diff --git a/MeshLib/MeshQuality/RadiusEdgeRatioMetric.cpp b/MeshLib/MeshQuality/RadiusEdgeRatioMetric.cpp index 4bf3257f3fefeb425f511164d710b2186a29ba90..97edcf556754a03700e13e7fe0837e4665e1f624 100644 --- a/MeshLib/MeshQuality/RadiusEdgeRatioMetric.cpp +++ b/MeshLib/MeshQuality/RadiusEdgeRatioMetric.cpp @@ -19,19 +19,19 @@ namespace MeshLib { - RadiusEdgeRatioMetric::RadiusEdgeRatioMetric(Mesh const& mesh) -: ElementQualityMetric(mesh) -{} + : ElementQualityMetric(mesh) +{ +} -void RadiusEdgeRatioMetric::calculateQuality () +void RadiusEdgeRatioMetric::calculateQuality() { std::vector<MeshLib::Element*> const& elements(_mesh.getElements()); - std::size_t const nElements (_mesh.getNumberOfElements()); + std::size_t const nElements(_mesh.getNumberOfElements()); for (std::size_t k(0); k < nElements; k++) { - Element const& elem (*elements[k]); - std::size_t const n_nodes (elem.getNumberOfBaseNodes()); + Element const& elem(*elements[k]); + std::size_t const n_nodes(elem.getNumberOfBaseNodes()); std::vector<MathLib::Point3d*> pnts(n_nodes); std::copy_n(elem.getNodes(), n_nodes, pnts.begin()); GeoLib::MinimalBoundingSphere const s(pnts); @@ -40,4 +40,4 @@ void RadiusEdgeRatioMetric::calculateQuality () } } -} // end namespace MeshLib +} // end namespace MeshLib diff --git a/MeshLib/MeshQuality/SizeDifferenceMetric.cpp b/MeshLib/MeshQuality/SizeDifferenceMetric.cpp index 645a7d7e0cbf5f2ca93c93f943008542d3937d6d..84965d1714f9d69193245390888596dbb0932fc9 100644 --- a/MeshLib/MeshQuality/SizeDifferenceMetric.cpp +++ b/MeshLib/MeshQuality/SizeDifferenceMetric.cpp @@ -18,39 +18,40 @@ namespace MeshLib { - -SizeDifferenceMetric::SizeDifferenceMetric(Mesh const& mesh) : -ElementQualityMetric(mesh) -{ } +SizeDifferenceMetric::SizeDifferenceMetric(Mesh const& mesh) + : ElementQualityMetric(mesh) +{ +} void SizeDifferenceMetric::calculateQuality() { std::vector<MeshLib::Element*> const& elements(_mesh.getElements()); - std::size_t const nElements (_mesh.getNumberOfElements()); - std::size_t const mesh_dim (_mesh.getDimension()); + std::size_t const nElements(_mesh.getNumberOfElements()); + std::size_t const mesh_dim(_mesh.getDimension()); - for (std::size_t k=0; k < nElements; ++k) + for (std::size_t k = 0; k < nElements; ++k) { - Element const& elem (*elements[k]); + Element const& elem(*elements[k]); if (elem.getDimension() < mesh_dim) { _element_quality_metric[k] = 0; continue; } - std::size_t const n_neighbors (elem.getNumberOfNeighbors()); - double const vol_a (elem.getContent()); + std::size_t const n_neighbors(elem.getNumberOfNeighbors()); + double const vol_a(elem.getContent()); double worst_ratio(1.0); - for (std::size_t i=0; i < n_neighbors; ++i) + for (std::size_t i = 0; i < n_neighbors; ++i) { - MeshLib::Element const*const neighbor (elem.getNeighbor(i)); + MeshLib::Element const* const neighbor(elem.getNeighbor(i)); if (neighbor == nullptr) { continue; } - double const vol_b (neighbor->getContent()); - double const ratio = (vol_a > vol_b) ? vol_b / vol_a : vol_a / vol_b; + double const vol_b(neighbor->getContent()); + double const ratio = + (vol_a > vol_b) ? vol_b / vol_a : vol_a / vol_b; if (ratio < worst_ratio) { worst_ratio = ratio; @@ -60,4 +61,4 @@ void SizeDifferenceMetric::calculateQuality() } } -} // end namespace MeshLib +} // end namespace MeshLib diff --git a/MeshLib/MeshSearch/ElementSearch.cpp b/MeshLib/MeshSearch/ElementSearch.cpp index 025460c8c58eb7f3f512fb6612a50a40b56e8528..7e5c8874790392bf3db578605f7347b03501cdfe 100644 --- a/MeshLib/MeshSearch/ElementSearch.cpp +++ b/MeshLib/MeshSearch/ElementSearch.cpp @@ -14,19 +14,17 @@ #include "MeshLib/Elements/Element.h" #include "MeshLib/Node.h" -namespace MeshLib { - -ElementSearch::ElementSearch(const MeshLib::Mesh &mesh) - : _mesh(mesh) +namespace MeshLib { -} +ElementSearch::ElementSearch(const MeshLib::Mesh& mesh) : _mesh(mesh) {} template <typename Container, typename Predicate> std::vector<std::size_t> filter(Container const& container, Predicate const& p) { std::vector<std::size_t> matchedIDs; std::size_t i = 0; - for (auto value : container) { + for (auto value : container) + { if (p(value)) { matchedIDs.push_back(i); @@ -38,8 +36,9 @@ std::vector<std::size_t> filter(Container const& container, Predicate const& p) std::size_t ElementSearch::searchByElementType(MeshElemType eleType) { - auto matchedIDs = filter(_mesh.getElements(), - [&](MeshLib::Element* e) { return e->getGeomType()==eleType; }); + auto matchedIDs = filter(_mesh.getElements(), [&](MeshLib::Element* e) { + return e->getGeomType() == eleType; + }); this->updateUnion(matchedIDs); return matchedIDs.size(); @@ -47,34 +46,34 @@ std::size_t ElementSearch::searchByElementType(MeshElemType eleType) std::size_t ElementSearch::searchByContent(double eps) { - auto matchedIDs = filter(_mesh.getElements(), - [&eps](MeshLib::Element* e) { return e->getContent() < eps; }); + auto matchedIDs = filter(_mesh.getElements(), [&eps](MeshLib::Element* e) { + return e->getContent() < eps; + }); this->updateUnion(matchedIDs); return matchedIDs.size(); } -std::size_t ElementSearch::searchByBoundingBox( - GeoLib::AABB const& aabb) +std::size_t ElementSearch::searchByBoundingBox(GeoLib::AABB const& aabb) { - auto matchedIDs = filter(_mesh.getElements(), - [&aabb](MeshLib::Element* e) { - std::size_t const nElemNodes (e->getNumberOfBaseNodes()); - for (std::size_t n = 0; n < nElemNodes; ++n) + auto matchedIDs = filter(_mesh.getElements(), [&aabb](MeshLib::Element* e) { + std::size_t const nElemNodes(e->getNumberOfBaseNodes()); + for (std::size_t n = 0; n < nElemNodes; ++n) + { + if (aabb.containsPoint(*e->getNode(n), 0)) { - if (aabb.containsPoint(*e->getNode(n), 0)) - { - return true; // any node of element is in aabb. - } + return true; // any node of element is in aabb. } - return false; // no nodes of element are in aabb. - }); + } + return false; // no nodes of element are in aabb. + }); this->updateUnion(matchedIDs); return matchedIDs.size(); } -std::size_t ElementSearch::searchByNodeIDs(const std::vector<std::size_t> &nodes) +std::size_t ElementSearch::searchByNodeIDs( + const std::vector<std::size_t>& nodes) { std::vector<std::size_t> connected_elements; for (std::size_t node_id : nodes) @@ -91,12 +90,13 @@ std::size_t ElementSearch::searchByNodeIDs(const std::vector<std::size_t> &nodes return connected_elements.size(); } -void ElementSearch::updateUnion(const std::vector<std::size_t> &vec) +void ElementSearch::updateUnion(const std::vector<std::size_t>& vec) { std::vector<std::size_t> vec_temp(vec.size() + _marked_elements.size()); - auto it = std::set_union(vec.begin(), vec.end(), _marked_elements.begin(), _marked_elements.end(), vec_temp.begin()); + auto it = std::set_union(vec.begin(), vec.end(), _marked_elements.begin(), + _marked_elements.end(), vec_temp.begin()); vec_temp.resize(it - vec_temp.begin()); _marked_elements.assign(vec_temp.begin(), vec_temp.end()); } -} // end namespace MeshLib +} // end namespace MeshLib diff --git a/MeshLib/MeshSearch/MeshElementGrid.cpp b/MeshLib/MeshSearch/MeshElementGrid.cpp index 0aa33c7554d144958547a50a731861ce8151b4d1..848cba2fc97bc2ec87a32883390dedd87f74e3b5 100644 --- a/MeshLib/MeshSearch/MeshElementGrid.cpp +++ b/MeshLib/MeshSearch/MeshElementGrid.cpp @@ -15,10 +15,9 @@ #include <cmath> #include <memory> +#include "../Elements/Element.h" #include "../Mesh.h" #include "../Node.h" -#include "../Elements/Element.h" - #include "GeoLib/GEOObjects.h" namespace MeshLib @@ -27,13 +26,14 @@ MeshElementGrid::MeshElementGrid(MeshLib::Mesh const& mesh) : _aabb{mesh.getNodes().cbegin(), mesh.getNodes().cend()}, _n_steps({{1, 1, 1}}) { - auto getDimensions = - [](MathLib::Point3d const& min, MathLib::Point3d const& max) - { + auto getDimensions = [](MathLib::Point3d const& min, + MathLib::Point3d const& max) { std::bitset<3> dim; // all bits are set to zero. - for (std::size_t k(0); k < 3; ++k) { + for (std::size_t k(0); k < 3; ++k) + { double const tolerance( - std::nexttoward(max[k],std::numeric_limits<double>::max())-max[k]); + std::nexttoward(max[k], std::numeric_limits<double>::max()) - + max[k]); if (std::abs(max[k] - min[k]) > tolerance) { dim[k] = true; @@ -46,58 +46,76 @@ MeshElementGrid::MeshElementGrid(MeshLib::Mesh const& mesh) MathLib::Point3d const& max_pnt(_aabb.getMaxPoint()); auto const dim = getDimensions(min_pnt, max_pnt); - std::array<double, 3> delta{{ max_pnt[0] - min_pnt[0], - max_pnt[1] - min_pnt[1], max_pnt[2] - min_pnt[2] }}; + std::array<double, 3> delta{{max_pnt[0] - min_pnt[0], + max_pnt[1] - min_pnt[1], + max_pnt[2] - min_pnt[2]}}; const std::size_t n_eles(mesh.getNumberOfElements()); const std::size_t n_eles_per_cell(100); // *** condition: n_eles / n_cells < n_eles_per_cell // where n_cells = _n_steps[0] * _n_steps[1] * _n_steps[2] - // *** with _n_steps[0] = ceil(pow(n_eles*delta[0]*delta[0]/(n_eles_per_cell*delta[1]*delta[2]), 1/3.))); + // *** with _n_steps[0] = + // ceil(pow(n_eles*delta[0]*delta[0]/(n_eles_per_cell*delta[1]*delta[2]), + // 1/3.))); // _n_steps[1] = _n_steps[0] * delta[1]/delta[0], // _n_steps[2] = _n_steps[0] * delta[2]/delta[0] - auto sc_ceil = [](double v){ + auto sc_ceil = [](double v) { return static_cast<std::size_t>(std::ceil(v)); }; - switch (dim.count()) { - case 3: // 3d case - _n_steps[0] = sc_ceil(std::cbrt( - n_eles*delta[0]*delta[0]/(n_eles_per_cell*delta[1]*delta[2]))); - _n_steps[1] = sc_ceil(_n_steps[0] * std::min(delta[1] / delta[0], 100.0)); - _n_steps[2] = sc_ceil(_n_steps[0] * std::min(delta[2] / delta[0], 100.0)); - break; - case 2: // 2d cases - if (dim[0] && dim[2]) { // 2d case: xz plane, y = const - _n_steps[0] = sc_ceil(std::sqrt(n_eles*delta[0]/(n_eles_per_cell*delta[2]))); - _n_steps[2] = sc_ceil(_n_steps[0]*delta[2]/delta[0]); - } - else if (dim[0] && dim[1]) { // 2d case: xy plane, z = const - _n_steps[0] = sc_ceil(std::sqrt(n_eles*delta[0]/(n_eles_per_cell*delta[1]))); - _n_steps[1] = sc_ceil(_n_steps[0] * delta[1]/delta[0]); - } - else if (dim[1] && dim[2]) { // 2d case: yz plane, x = const - _n_steps[1] = sc_ceil(std::sqrt(n_eles*delta[1]/(n_eles_per_cell*delta[2]))); - _n_steps[2] = sc_ceil(n_eles * delta[2] / (n_eles_per_cell*delta[1])); - } - break; - case 1: // 1d cases - for (std::size_t k(0); k<3; ++k) { - if (dim[k]) { - _n_steps[k] = sc_ceil(static_cast<double>(n_eles)/n_eles_per_cell); + switch (dim.count()) + { + case 3: // 3d case + _n_steps[0] = + sc_ceil(std::cbrt(n_eles * delta[0] * delta[0] / + (n_eles_per_cell * delta[1] * delta[2]))); + _n_steps[1] = + sc_ceil(_n_steps[0] * std::min(delta[1] / delta[0], 100.0)); + _n_steps[2] = + sc_ceil(_n_steps[0] * std::min(delta[2] / delta[0], 100.0)); + break; + case 2: // 2d cases + if (dim[0] && dim[2]) + { // 2d case: xz plane, y = const + _n_steps[0] = sc_ceil(std::sqrt(n_eles * delta[0] / + (n_eles_per_cell * delta[2]))); + _n_steps[2] = sc_ceil(_n_steps[0] * delta[2] / delta[0]); + } + else if (dim[0] && dim[1]) + { // 2d case: xy plane, z = const + _n_steps[0] = sc_ceil(std::sqrt(n_eles * delta[0] / + (n_eles_per_cell * delta[1]))); + _n_steps[1] = sc_ceil(_n_steps[0] * delta[1] / delta[0]); + } + else if (dim[1] && dim[2]) + { // 2d case: yz plane, x = const + _n_steps[1] = sc_ceil(std::sqrt(n_eles * delta[1] / + (n_eles_per_cell * delta[2]))); + _n_steps[2] = + sc_ceil(n_eles * delta[2] / (n_eles_per_cell * delta[1])); + } + break; + case 1: // 1d cases + for (std::size_t k(0); k < 3; ++k) + { + if (dim[k]) + { + _n_steps[k] = + sc_ceil(static_cast<double>(n_eles) / n_eles_per_cell); + } } - } } // some frequently used expressions to fill the vector of elements per grid // cell - for (std::size_t k(0); k<3; k++) { + for (std::size_t k(0); k < 3; k++) + { _step_sizes[k] = delta[k] / _n_steps[k]; _inverse_step_sizes[k] = 1.0 / _step_sizes[k]; } - _elements_in_grid_box.resize(_n_steps[0]*_n_steps[1]*_n_steps[2]); + _elements_in_grid_box.resize(_n_steps[0] * _n_steps[1] * _n_steps[2]); sortElementsInGridCells(mesh); } @@ -115,7 +133,8 @@ void MeshElementGrid::sortElementsInGridCells(MeshLib::Mesh const& mesh) { for (auto const element : mesh.getElements()) { - if (! sortElementInGridCells(*element)) { + if (!sortElementInGridCells(*element)) + { OGS_FATAL("Sorting element (id={:d}) into mesh element grid.", element->getID()); } @@ -124,20 +143,25 @@ void MeshElementGrid::sortElementsInGridCells(MeshLib::Mesh const& mesh) bool MeshElementGrid::sortElementInGridCells(MeshLib::Element const& element) { - std::array<std::size_t,3> min{}; - std::array<std::size_t,3> max{}; - std::pair<bool, std::array<std::size_t, 3>> c( - getGridCellCoordinates(*(static_cast<MathLib::Point3d const*>(element.getNode(0))))); - if (c.first) { + std::array<std::size_t, 3> min{}; + std::array<std::size_t, 3> max{}; + std::pair<bool, std::array<std::size_t, 3>> c(getGridCellCoordinates( + *(static_cast<MathLib::Point3d const*>(element.getNode(0))))); + if (c.first) + { min = c.second; max = min; - } else { + } + else + { return false; } - for (std::size_t k(1); k<element.getNumberOfNodes(); ++k) { + for (std::size_t k(1); k < element.getNumberOfNodes(); ++k) + { // compute coordinates of the grid for each node of the element - c = getGridCellCoordinates(*(static_cast<MathLib::Point3d const*>(element.getNode(k)))); + c = getGridCellCoordinates( + *(static_cast<MathLib::Point3d const*>(element.getNode(k)))); if (!c.first) { return false; @@ -156,7 +180,7 @@ bool MeshElementGrid::sortElementInGridCells(MeshLib::Element const& element) } } - const std::size_t n_plane(_n_steps[0]*_n_steps[1]); + const std::size_t n_plane(_n_steps[0] * _n_steps[1]); // If a node of an element is almost equal to the upper right point of the // AABB the grid cell coordinates computed by getGridCellCoordintes() could @@ -168,10 +192,14 @@ bool MeshElementGrid::sortElementInGridCells(MeshLib::Element const& element) } // insert the element into the grid cells - for (std::size_t i(min[0]); i<=max[0]; i++) { - for (std::size_t j(min[1]); j<=max[1]; j++) { - for (std::size_t k(min[2]); k<=max[2]; k++) { - _elements_in_grid_box[i+j*_n_steps[0]+k*n_plane].push_back(&element); + for (std::size_t i(min[0]); i <= max[0]; i++) + { + for (std::size_t j(min[1]); j <= max[1]; j++) + { + for (std::size_t k(min[2]); k <= max[2]; k++) + { + _elements_in_grid_box[i + j * _n_steps[0] + k * n_plane] + .push_back(&element); } } } @@ -185,15 +213,21 @@ MeshElementGrid::getGridCellCoordinates(MathLib::Point3d const& p) const bool valid(true); std::array<std::size_t, 3> coords{}; - for (std::size_t k(0); k<3; ++k) { - const double d(p[k]-_aabb.getMinPoint()[k]); - if (d < 0.0) { + for (std::size_t k(0); k < 3; ++k) + { + const double d(p[k] - _aabb.getMinPoint()[k]); + if (d < 0.0) + { valid = false; coords[k] = 0; - } else if (_aabb.getMaxPoint()[k] <= p[k]) { + } + else if (_aabb.getMaxPoint()[k] <= p[k]) + { valid = false; - coords[k] = _n_steps[k]-1; - } else { + coords[k] = _n_steps[k] - 1; + } + else + { coords[k] = static_cast<std::size_t>(d * _inverse_step_sizes[k]); } } diff --git a/MeshLib/MeshSearch/NodeSearch.cpp b/MeshLib/MeshSearch/NodeSearch.cpp index 68f5034f75250921a16fed7fd3d431b6b5e8158e..d5d5b58c7db5f29150460cfefb0f6a3c242f0274 100644 --- a/MeshLib/MeshSearch/NodeSearch.cpp +++ b/MeshLib/MeshSearch/NodeSearch.cpp @@ -13,19 +13,16 @@ #include <memory> #include <set> +#include "MeshLib/Elements/Element.h" #include "MeshLib/Mesh.h" #include "MeshLib/Node.h" -#include "MeshLib/Elements/Element.h" -namespace MeshLib { - -NodeSearch::NodeSearch(const MeshLib::Mesh &mesh) - : _mesh(mesh) +namespace MeshLib { -} +NodeSearch::NodeSearch(const MeshLib::Mesh& mesh) : _mesh(mesh) {} std::size_t NodeSearch::searchNodesConnectedToOnlyGivenElements( - const std::vector<std::size_t> &elements) + const std::vector<std::size_t>& elements) { // Find out by how many elements a node would be removed. // @@ -33,19 +30,19 @@ std::size_t NodeSearch::searchNodesConnectedToOnlyGivenElements( // algorithm might be more memory efficient. std::vector<std::size_t> node_marked_counts(_mesh.getNumberOfNodes(), 0); - for(std::size_t eid : elements) + for (std::size_t eid : elements) { auto* e = _mesh.getElement(eid); - for (unsigned i=0; i<e->getNumberOfNodes(); i++) { + for (unsigned i = 0; i < e->getNumberOfNodes(); i++) + { node_marked_counts[e->getNodeIndex(i)]++; } } - // Push back nodes which counts are equal to number of connected elements to // that node. std::vector<std::size_t> connected_nodes; - for (std::size_t i=0; i<node_marked_counts.size(); i++) + for (std::size_t i = 0; i < node_marked_counts.size(); i++) { if (node_marked_counts[i] == _mesh.getNode(i)->getElements().size()) { @@ -59,8 +56,8 @@ std::size_t NodeSearch::searchNodesConnectedToOnlyGivenElements( std::size_t NodeSearch::searchUnused() { - const std::size_t nNodes (_mesh.getNumberOfNodes()); - const std::vector<MeshLib::Node*> &nodes (_mesh.getNodes()); + const std::size_t nNodes(_mesh.getNumberOfNodes()); + const std::vector<MeshLib::Node*>& nodes(_mesh.getNodes()); std::vector<std::size_t> del_node_idx; for (unsigned i = 0; i < nNodes; ++i) @@ -101,8 +98,8 @@ std::size_t NodeSearch::searchBoundaryNodes() continue; } - std::size_t const n_edges (elem->getNumberOfEdges()); - for (std::size_t i=0; i<n_edges; ++i) + std::size_t const n_edges(elem->getNumberOfEdges()); + for (std::size_t i = 0; i < n_edges; ++i) { if (elem->getNeighbor(i) != nullptr) { @@ -129,8 +126,8 @@ std::size_t NodeSearch::searchBoundaryNodes() continue; } - std::size_t const n_faces (elem->getNumberOfFaces()); - for (std::size_t i=0; i<n_faces; ++i) + std::size_t const n_faces(elem->getNumberOfFaces()); + for (std::size_t i = 0; i < n_faces; ++i) { if (elem->getNeighbor(i) != nullptr) { @@ -145,17 +142,19 @@ std::size_t NodeSearch::searchBoundaryNodes() } } std::sort(vec_boundary_nodes.begin(), vec_boundary_nodes.end()); - vec_boundary_nodes.erase(std::unique(vec_boundary_nodes.begin(), vec_boundary_nodes.end()), vec_boundary_nodes.end()); - + vec_boundary_nodes.erase( + std::unique(vec_boundary_nodes.begin(), vec_boundary_nodes.end()), + vec_boundary_nodes.end()); this->updateUnion(vec_boundary_nodes); return vec_boundary_nodes.size(); } -void NodeSearch::updateUnion(const std::vector<std::size_t> &vec) +void NodeSearch::updateUnion(const std::vector<std::size_t>& vec) { std::vector<std::size_t> vec_temp(vec.size() + _marked_nodes.size()); - auto it = std::set_union(vec.begin(), vec.end(), _marked_nodes.begin(), _marked_nodes.end(), vec_temp.begin()); + auto it = std::set_union(vec.begin(), vec.end(), _marked_nodes.begin(), + _marked_nodes.end(), vec_temp.begin()); vec_temp.resize(it - vec_temp.begin()); _marked_nodes.assign(vec_temp.begin(), vec_temp.end()); } @@ -173,10 +172,9 @@ std::vector<Node*> getUniqueNodes(std::vector<Element*> const& elements) std::vector<Node*> nodes; nodes.reserve(nodes_set.size()); - std::move(nodes_set.cbegin(), nodes_set.cend(), - std::back_inserter(nodes)); + std::move(nodes_set.cbegin(), nodes_set.cend(), std::back_inserter(nodes)); return nodes; } -} // end namespace MeshLib +} // end namespace MeshLib diff --git a/MeshLib/MeshSurfaceExtraction.cpp b/MeshLib/MeshSurfaceExtraction.cpp index 13939a2f9743f219bda400d5dbc7877e5170f79c..3d37732aa82c7aa7f17014c91eb07e24f310073a 100644 --- a/MeshLib/MeshSurfaceExtraction.cpp +++ b/MeshLib/MeshSurfaceExtraction.cpp @@ -15,9 +15,9 @@ #include "MeshSurfaceExtraction.h" #include <boost/math/constants/constants.hpp> -#include "BaseLib/Logging.h" #include <memory> +#include "BaseLib/Logging.h" #include "MeshLib/Elements/Line.h" #include "MeshLib/Elements/Point.h" #include "MeshLib/Elements/Quad.h" diff --git a/MeshLib/Node.cpp b/MeshLib/Node.cpp index 886d16de7b7874c2cdbc249edd32d610bad07390..df1d730dde44220affa41c4bbabf6a53b0cf930a 100644 --- a/MeshLib/Node.cpp +++ b/MeshLib/Node.cpp @@ -13,13 +13,14 @@ */ #include "MeshLib/Node.h" -#include "Elements/Element.h" -namespace MeshLib { +#include "Elements/Element.h" +namespace MeshLib +{ Node::Node(const double coords[3], std::size_t id) : MathLib::Point3dWithID( - std::array<double,3>{{coords[0], coords[1], coords[2]}}, id) + std::array<double, 3>{{coords[0], coords[1], coords[2]}}, id) { } @@ -29,14 +30,11 @@ Node::Node(std::array<double, 3> const& coords, std::size_t id) } Node::Node(double x, double y, double z, std::size_t id) - : MathLib::Point3dWithID(std::array<double,3>({{x, y, z}}), id) + : MathLib::Point3dWithID(std::array<double, 3>({{x, y, z}}), id) { } -Node::Node(const Node &node) - : MathLib::Point3dWithID(node._x, node.getID()) -{ -} +Node::Node(const Node& node) : MathLib::Point3dWithID(node._x, node.getID()) {} void Node::updateCoordinates(double x, double y, double z) { @@ -44,7 +42,7 @@ void Node::updateCoordinates(double x, double y, double z) _x[1] = y; _x[2] = z; - const std::size_t nElements (this->_elements.size()); + const std::size_t nElements(this->_elements.size()); for (std::size_t i = 0; i < nElements; i++) { _elements[i]->computeVolume(); diff --git a/MeshLib/Properties.cpp b/MeshLib/Properties.cpp index ef5ae723f2a8b319fa8339df485be05edadc4013..1def95c05236d0fcdbf4f8db8d2eb68ed57ce9c2 100644 --- a/MeshLib/Properties.cpp +++ b/MeshLib/Properties.cpp @@ -14,13 +14,12 @@ namespace MeshLib { - void Properties::removePropertyVector(std::string const& name) { std::map<std::string, PropertyVectorBase*>::const_iterator it( - _properties.find(name) - ); - if (it == _properties.end()) { + _properties.find(name)); + if (it == _properties.end()) + { WARN("A property of the name '{:s}' does not exist.", name); return; } diff --git a/MeshLib/Vtk/VtkMappedMeshSource.cpp b/MeshLib/Vtk/VtkMappedMeshSource.cpp index b24b321cd45322cc97cdfa12a51b88a315451f3a..f99251321bdf855dbb2cd2768f246271d7d33c42 100644 --- a/MeshLib/Vtk/VtkMappedMeshSource.cpp +++ b/MeshLib/Vtk/VtkMappedMeshSource.cpp @@ -13,8 +13,6 @@ */ #include "VtkMappedMeshSource.h" -#include <vector> - #include <vtkCellType.h> #include <vtkDemandDrivenPipeline.h> #include <vtkInformation.h> @@ -22,6 +20,8 @@ #include <vtkSmartPointer.h> #include <vtkStreamingDemandDrivenPipeline.h> +#include <vector> + #include "MeshLib/Elements/Element.h" #include "MeshLib/Mesh.h" #include "MeshLib/VtkOGSEnum.h" @@ -172,8 +172,8 @@ int VtkMappedMeshSource::RequestData(vtkInformation* /*request*/, { addProperty(*p); } - else if (auto p = - dynamic_cast<PropertyVector<unsigned long long>*>(property)) + else if (auto p = dynamic_cast<PropertyVector<unsigned long long>*>( + property)) { addProperty(*p); } diff --git a/MeshLib/VtkOGSEnum.cpp b/MeshLib/VtkOGSEnum.cpp index 55ed23e1cec7b8fe4d729dcd211d7dd5fdcfd3d9..55d91c3c340774c02482ef81d11fe239cd9c596e 100644 --- a/MeshLib/VtkOGSEnum.cpp +++ b/MeshLib/VtkOGSEnum.cpp @@ -10,10 +10,10 @@ #include "VtkOGSEnum.h" -#include "BaseLib/Error.h" - #include <vtkCellType.h> +#include "BaseLib/Error.h" + int OGSToVtkCellType(MeshLib::CellType ogs) { switch (ogs) @@ -63,5 +63,3 @@ int OGSToVtkCellType(MeshLib::CellType ogs) ogs); } } - - diff --git a/MeshLib/convertMeshToGeo.cpp b/MeshLib/convertMeshToGeo.cpp index 65ff3b149bcc91229600b26f96f42281df93345a..c3fde08fffacbd81bb25f6b266e60d7ac1f79e90 100644 --- a/MeshLib/convertMeshToGeo.cpp +++ b/MeshLib/convertMeshToGeo.cpp @@ -15,13 +15,11 @@ #include "convertMeshToGeo.h" #include "BaseLib/Logging.h" - +#include "Elements/Quad.h" +#include "Elements/Tri.h" #include "GeoLib/GEOObjects.h" #include "GeoLib/Surface.h" #include "GeoLib/Triangle.h" - -#include "Elements/Quad.h" -#include "Elements/Tri.h" #include "Mesh.h" #include "MeshEditing/MeshRevision.h" #include "MeshInformation.h" diff --git a/MeshLib/findElementsWithinRadius.cpp b/MeshLib/findElementsWithinRadius.cpp index 36dd45ef9ff0703fcbb4c460348b0a346692d442..e40a104524d1acd36f208c0922758b4ff056152e 100644 --- a/MeshLib/findElementsWithinRadius.cpp +++ b/MeshLib/findElementsWithinRadius.cpp @@ -14,12 +14,12 @@ #include <unordered_set> #include <vector> -#include "MathLib/MathTools.h" - #include "Elements/Element.h" +#include "MathLib/MathTools.h" #include "Node.h" -namespace MeshLib { +namespace MeshLib +{ std::vector<std::size_t> findElementsWithinRadius(Element const& start_element, double const radius_squared) {