diff --git a/FileIO/MeshIO.cpp b/FileIO/MeshIO.cpp index 7ea88cbbdebf4ea2add0bfb0bcc1eefc38892892..6fe6ffae4081afb84d647c25ef0eeb8a063e03b8 100644 --- a/FileIO/MeshIO.cpp +++ b/FileIO/MeshIO.cpp @@ -173,10 +173,7 @@ MeshLib::Element* MeshIO::readElement(const std::string& line, const std::vector int MeshIO::write(std::ostream &out) { - (void)out; - /* - if(!_mesh) - { + if(!_mesh) { std::cout << "OGSMeshIO cannot write: no mesh set!" << std::endl; return 0; } @@ -185,25 +182,22 @@ int MeshIO::write(std::ostream &out) out << "#FEM_MSH" << std::endl; - out << "$PCS_TYPE" << std::endl << " " << _mesh->pcs_name << std::endl; + out << "$PCS_TYPE" << std::endl << " NO_PCS" << std::endl; out << "$NODES" << std::endl << " "; - const size_t n_nodes(_mesh->GetNodesNumber(false)); + const size_t n_nodes(_mesh->getNNodes()); out << n_nodes << std::endl; - for (size_t i(0); i < n_nodes; i++) - { - double const* const coords (_mesh->nod_vector[i]->getCoords()); - out << i << " " << coords[0] << " " << coords[1] << " " << coords[2] << std::endl; + for (size_t i(0); i < n_nodes; i++) { + out << i << " " << *(_mesh->getNode(i)) << std::endl; } out << "$ELEMENTS" << std::endl << " "; - writeElementsExceptLines(_mesh->ele_vector, out); + writeElementsExceptLines(_mesh->getElements(), out); out << " $LAYER" << std::endl; - out << " "; - out << _mesh->_n_msh_layer << std::endl; + out << " 0" << std::endl; out << "#STOP" << std::endl; - */ + return 1; } @@ -212,4 +206,37 @@ void MeshIO::setMesh(const MeshLib::Mesh* mesh) _mesh = mesh; } +void MeshIO::writeElementsExceptLines(std::vector<MeshLib::Element*> const& ele_vec, std::ostream &out) +{ + const size_t ele_vector_size (ele_vec.size()); + const double epsilon (std::numeric_limits<double>::epsilon()); + std::vector<bool> non_line_element (ele_vector_size, true); + std::vector<bool> non_null_element (ele_vector_size, true); + size_t n_elements(0); + + for (size_t i(0); i < ele_vector_size; i++) { + if ((ele_vec[i])->getType() == MshElemType::LINE) { + non_line_element[i] = false; + non_null_element[i] = false; + } else { + if (ele_vec[i]->getContent() < epsilon) { + non_null_element[i] = false; + } else { + n_elements++; + } + } + } + out << n_elements << std::endl; + for (size_t i(0), k(0); i < ele_vector_size; i++) { + if (non_line_element[i] && non_null_element[i]) { + out << k << " 0 " << MshElemType2String(ele_vec[i]->getType()) << " "; + for(size_t j = 0; j < ele_vec[i]->getNNodes()-1; j++) + out << ele_vec[i]->getNode(j)->getID() << " "; + out << ele_vec[i]->getNode(ele_vec[i]->getNNodes()-1)->getID() << std::endl; + k++; + } + } +} + + } // end namespace FileIO diff --git a/FileIO/MeshIO.h b/FileIO/MeshIO.h index 6472e4fdefca01a1919a2aa2e142b766bd27e463..d78462b0bd3f23cc9ee406f027f7cddd660291b2 100644 --- a/FileIO/MeshIO.h +++ b/FileIO/MeshIO.h @@ -51,6 +51,7 @@ protected: int write(std::ostream &out); private: + void writeElementsExceptLines(std::vector<MeshLib::Element*> const& ele_vec, std::ostream &out); MeshLib::Element* readElement(const std::string& line, const std::vector<MeshLib::Node*> &nodes); double* _edge_length[2]; diff --git a/GeoLib/Grid.h b/GeoLib/Grid.h index acb1b5b51b6838caf2edb558972f625aa3ec9f42..b5c5f0f7cbeaf991feac6da4cdbcec50fd92b4db 100644 --- a/GeoLib/Grid.h +++ b/GeoLib/Grid.h @@ -17,6 +17,7 @@ // GeoLib #include "AxisAlignedBoundingBox.h" +#include "GEOObjects.h" #ifndef NDEBUG // BaseLib @@ -30,8 +31,8 @@ class Grid : public GeoLib::AABB { public: /** - * @brief The constructor of the grid object takes a vector of points. Furthermore the - * use can specify the maximum number of points per grid cell (in the average!). + * @brief The constructor of the grid object takes a vector of points or nodes. Furthermore the + * use can specify the *average* maximum number of points per grid cell. * * The number of grid cells are computed with the following formula * \f$\frac{n_{points}}{n_{cells}} \le n_{max\_per\_cell}\f$ @@ -54,9 +55,9 @@ public: } /** - * The method calculates the grid cell within the grid to given point is belonging to, i.e., - * the (internal) coordinates of the grid are computed. The method searches the actual - * grid cell and all its possible neighbors for the POINT object which has the smallest + * The method calculates the grid cell the given point is belonging to, i.e., + * the (internal) coordinates of the grid cell are computed. The method searches the actual + * grid cell and all its neighbors for the POINT object which has the smallest * distance. A pointer to this object is returned. * * If there is not such a point, i.e., all the searched grid cells do not contain any @@ -68,15 +69,15 @@ public: POINT const* getNearestPoint(double const*const pnt) const; /** - * Method fetchs all points that are located within grid cells that intersects - * the axis aligned cube defined by its center and half edge length. + * Method fetches the vectors of all grid cells intersecting the axis aligned cube + * defined by its center and half edge length. * * @param pnt (input) the center point of the axis aligned cube * @param half_len (input) half of the edge length of the axis aligned cube - * @param pnts (output) all points within grid cells that intersects + * @param pnts (output) vector of vectors of points within grid cells that intersects * the axis aligned cube */ - void getPointsWithinCube(double const*const pnt, double half_len, std::vector<POINT*>& pnts) const; + void getVecsOfGridCellsIntersectingCube(double const*const pnt, double half_len, std::vector<std::vector<POINT*> const*>& pnts) const; #ifndef NDEBUG /** @@ -289,15 +290,19 @@ POINT const* Grid<POINT>::getNearestPoint(double const*const pnt) const double len (sqrt(MathLib::sqrDist(pnt, nearest_pnt->getCoords()))); // search all other grid cells within the cube with the edge nodes - std::vector<POINT*> pnts; - getPointsWithinCube(pnt, len, pnts); - - const size_t n_pnts(pnts.size()); - for (size_t k(0); k<n_pnts; k++) { - const double sqr_dist (MathLib::sqrDist(pnt, pnts[k]->getCoords())); - if (sqr_dist < sqr_min_dist) { - sqr_min_dist = sqr_dist; - nearest_pnt = pnts[k]; + std::vector<std::vector<POINT*> const*> vecs_of_pnts; + getVecsOfGridCellsIntersectingCube(pnt, len, vecs_of_pnts); + + const size_t n_vecs(vecs_of_pnts.size()); + for (size_t j(0); j<n_vecs; j++) { + std::vector<POINT*> const& pnts(vecs_of_pnts[j]); + const size_t n_pnts(pnts.size()); + for (size_t k(0); k<n_pnts; k++) { + const double sqr_dist (MathLib::sqrDist(pnt, pnts[k]->getCoords())); + if (sqr_dist < sqr_min_dist) { + sqr_min_dist = sqr_dist; + nearest_pnt = pnts[k]; + } } } @@ -305,7 +310,7 @@ POINT const* Grid<POINT>::getNearestPoint(double const*const pnt) const } template <typename POINT> -void Grid<POINT>::getPointsWithinCube(double const*const pnt, double half_len, std::vector<POINT*>& pnts) const +void Grid<POINT>::getVecsOfGridCellsIntersectingCube(double const*const pnt, double half_len, std::vector<std::vector<POINT*> const*>& pnts) const { double tmp_pnt[3] = {pnt[0]-half_len, pnt[1]-half_len, pnt[2]-half_len}; // min size_t min_coords[3]; @@ -322,11 +327,12 @@ void Grid<POINT>::getPointsWithinCube(double const*const pnt, double half_len, s for (coords[1] = min_coords[1]; coords[1] < max_coords[1]+1; coords[1]++) { const size_t coords0_p_coords1_x_steps0 (coords[0] + coords[1] * _n_steps[0]); for (coords[2] = min_coords[2]; coords[2] < max_coords[2]+1; coords[2]++) { - // copy pnts - const size_t n_pnts(_grid_quad_to_node_map[coords0_p_coords1_x_steps0 + coords[2] * steps0_x_steps1].size()); - for (size_t k(0); k<n_pnts; k++) { - pnts.push_back(_grid_quad_to_node_map[coords0_p_coords1_x_steps0 + coords[2] * steps0_x_steps1][k]); - } + pnts.push_back (&(_grid_quad_to_node_map[coords0_p_coords1_x_steps0 + coords[2] * steps0_x_steps1])); +// // copy pnts +// const size_t n_pnts(_grid_quad_to_node_map[coords0_p_coords1_x_steps0 + coords[2] * steps0_x_steps1].size()); +// for (size_t k(0); k<n_pnts; k++) { +// pnts.push_back(_grid_quad_to_node_map[coords0_p_coords1_x_steps0 + coords[2] * steps0_x_steps1][k]); +// } } } } diff --git a/MathLib/LinAlg/Sparse/NestedDissectionPermutation/Cluster.cpp b/MathLib/LinAlg/Sparse/NestedDissectionPermutation/Cluster.cpp index bd6c8c77f055d914bdd09b68da6b7488e4127395..45054a8f35b9e1dd9fe53fd26866b64fc5cb205c 100644 --- a/MathLib/LinAlg/Sparse/NestedDissectionPermutation/Cluster.cpp +++ b/MathLib/LinAlg/Sparse/NestedDissectionPermutation/Cluster.cpp @@ -5,7 +5,7 @@ * http://www.opengeosys.com/LICENSE.txt * * - * \fileCluster.cpp + * \file Cluster.cpp * * Created on 2012-01-02 by Thomas Fischer */ diff --git a/MeshLib/Elements/Cell.cpp b/MeshLib/Elements/Cell.cpp index 7c1a9f0f5881aa79602c793d5aa9709c90f84233..47a1e7a8d38b74910f1f739fd295f60b9866e020 100644 --- a/MeshLib/Elements/Cell.cpp +++ b/MeshLib/Elements/Cell.cpp @@ -19,15 +19,13 @@ Cell::Cell(Node** nodes, MshElemType::type type, unsigned value) { } */ -Cell::Cell(MshElemType::type type, unsigned value) - : Element(type, value) +Cell::Cell(unsigned value) + : Element(value) { } Cell::~Cell() -{ - delete[] this->_neighbors; -} +{} } diff --git a/MeshLib/Elements/Cell.h b/MeshLib/Elements/Cell.h index d3d02df3fb71514c16073853329da7b5783d2d4e..32d4ac88176f55270d781bcedf1657c5e2d6a143 100644 --- a/MeshLib/Elements/Cell.h +++ b/MeshLib/Elements/Cell.h @@ -35,13 +35,20 @@ public: /// Destructor virtual ~Cell(); + /** + * This method is pure virtual and is inherited from class @sa Element. + * It has to be implemented in the derived classes of class Cell! + * @return a copy of the object + */ + virtual Element* clone() const = 0; + protected: /* /// Constructor for a generic mesh element containing an array of mesh nodes. Cell(Node** nodes, MshElemType::type type, unsigned value = 0); */ /// Constructor for a generic mesh element without an array of mesh nodes. - Cell(MshElemType::type type, unsigned value = 0); + Cell(unsigned value = 0); /// Calculate the volume of this 3d element. virtual double computeVolume() = 0; diff --git a/MeshLib/Elements/Edge.cpp b/MeshLib/Elements/Edge.cpp index efbb774d2b08d22acde5e29da8b83e78e12d23c9..bef135b9c763c2d85091a6466653ad0fc1a73050 100644 --- a/MeshLib/Elements/Edge.cpp +++ b/MeshLib/Elements/Edge.cpp @@ -18,24 +18,23 @@ namespace MeshLib { Edge::Edge(Node* nodes[2], unsigned value) - : Element(MshElemType::LINE, value) + : Element(value) { _nodes = nodes; this->_length = this->computeLength(); } Edge::Edge(Node* n0, Node* n1, unsigned value) - : Element(MshElemType::LINE, value) + : Element(value) { _nodes = new Node*[2]; _nodes[0] = n0; _nodes[1] = n1; - this->_length = this->computeLength(); } Edge::Edge(const Edge &edge) - : Element(MshElemType::LINE, edge.getValue()) + : Element(edge.getValue()) { _nodes = new Node*[2]; _nodes[0] = edge._nodes[0]; @@ -52,5 +51,10 @@ double Edge::computeLength() return sqrt(MathLib::sqrDist(_nodes[0]->getCoords(), _nodes[1]->getCoords())); } +Element* Edge::clone() const +{ + return new Edge(*this); +} + } diff --git a/MeshLib/Elements/Edge.h b/MeshLib/Elements/Edge.h index 2863dee83537203c4c3764b8227605d5904ce6ef..eea222fa8a362dc6766c1d751f87bab2a80d10c7 100644 --- a/MeshLib/Elements/Edge.h +++ b/MeshLib/Elements/Edge.h @@ -76,6 +76,9 @@ public: /// Get the number of nodes for this element. virtual unsigned getNNodes() const { return 2; }; + virtual MshElemType::type getType() const { return MshElemType::LINE; } + + virtual Element* clone() const; protected: /// Calculate the length of this 1d element. diff --git a/MeshLib/Elements/Element.cpp b/MeshLib/Elements/Element.cpp index b952c3fe319d171604fd86596e69f4f7584a0072..5541d038162357e2a8ffdeb129a86fbc31c9c606 100644 --- a/MeshLib/Elements/Element.cpp +++ b/MeshLib/Elements/Element.cpp @@ -23,14 +23,15 @@ Element::Element(Node** nodes, MshElemType::type type, unsigned value) { } */ -Element::Element(MshElemType::type type, unsigned value) - : _nodes(NULL), _type(type), _value(value) +Element::Element(unsigned value) + : _nodes(NULL), _value(value), _neighbors(NULL) { } Element::~Element() { - delete[] this->_nodes; + delete [] this->_nodes; + delete [] this->_neighbors; } const Element* Element::getEdge(unsigned i) const @@ -71,10 +72,16 @@ const Node* Element::getNode(unsigned i) const { if (i < getNNodes()) return _nodes[i]; - std::cerr << "Error in MeshLib::Element::getNode() - Index does not exist." << std::endl; + std::cerr << "Error in MeshLib::Element::getNode() - Index " << i << " in " << MshElemType2String(getType()) << " does not exist." << std::endl; return NULL; } +void Element::setNode(unsigned idx, Node* node) +{ + if (idx < getNNodes()) + _nodes[idx] = node; +} + unsigned Element::getNodeIndex(unsigned i) const { if (i<getNNodes()) diff --git a/MeshLib/Elements/Element.h b/MeshLib/Elements/Element.h index 33abe67e144a72ee17049f180985a86715c9609e..e108525699ef2733c413e9139d0eb98e991f7d90 100644 --- a/MeshLib/Elements/Element.h +++ b/MeshLib/Elements/Element.h @@ -41,6 +41,13 @@ public: /// Get node with local index i. const Node* getNode(unsigned i) const; + /** + * (Re)Sets the node of the element. + * @param idx the index of the pointer to a node within the element + * @param node a pointer to a node + */ + void setNode(unsigned idx, Node* node); + /// Get array of element nodes. Node* const* getNodes() const { return _nodes; }; @@ -75,7 +82,7 @@ public: unsigned getNodeIndex(unsigned i) const; /// Get the type of the mesh element (as a MshElemType-enum). - MshElemType::type getType() const { return _type; }; + virtual MshElemType::type getType() const = 0; /// Get the value for this element. unsigned getValue() const { return _value; }; @@ -85,15 +92,21 @@ public: /// Destructor virtual ~Element(); + /** + * Method clone is a pure virtual method in the abstract base class Element. + * It has to be implemented in the derived classes (for instance in class Hex). + * @return an exact copy of the object + */ + virtual Element* clone() const = 0; + protected: /// Constructor for a generic mesh element without an array of mesh nodes. - Element(MshElemType::type type, unsigned value = 0); + Element(unsigned value = 0); /// Return a specific edge node. virtual Node* getEdgeNode(unsigned edge_id, unsigned node_id) const = 0; Node** _nodes; - MshElemType::type _type; unsigned _value; Element** _neighbors; diff --git a/MeshLib/Elements/Face.cpp b/MeshLib/Elements/Face.cpp index cab9b9ece368a11d3f7894335f0b425982e055d6..751e3807d5d28a1323c3b8473829a102a43db72e 100644 --- a/MeshLib/Elements/Face.cpp +++ b/MeshLib/Elements/Face.cpp @@ -20,15 +20,13 @@ Face::Face(Node** nodes, MshElemType::type type, unsigned value) { } */ -Face::Face(MshElemType::type type, unsigned value) - : Element(type, value) +Face::Face(unsigned value) + : Element(value) { } Face::~Face() -{ - delete[] this->_neighbors; -} +{} diff --git a/MeshLib/Elements/Face.h b/MeshLib/Elements/Face.h index f3866aebd07a0e09ea09da2adb5bad77e323d3da..1d6347a1e53c46425fededdf7652609162ef2a05 100644 --- a/MeshLib/Elements/Face.h +++ b/MeshLib/Elements/Face.h @@ -44,13 +44,20 @@ public: /// Destructor virtual ~Face(); + /** + * This method is pure virtual and is inherited from class @sa Element. + * It has to be implemented in the derived classes of class Face! + * @return a copy of the object + */ + virtual Element* clone() const = 0; + protected: /* /// Constructor for a generic mesh element containing an array of mesh nodes. Face(Node** nodes, MshElemType::type type, unsigned value = 0); */ /// Constructor for a generic mesh element without an array of mesh nodes. - Face(MshElemType::type type, unsigned value = 0); + Face(unsigned value = 0); /// Calculate the area of this 2d element. virtual double computeArea() = 0; diff --git a/MeshLib/Elements/Hex.cpp b/MeshLib/Elements/Hex.cpp index 5d219578a7fcdb83fd5bb05dd3b443df25bc9ec6..f296292fd41373c18448b7070fee3a2cf497b696 100644 --- a/MeshLib/Elements/Hex.cpp +++ b/MeshLib/Elements/Hex.cpp @@ -47,7 +47,7 @@ const unsigned Hex::_edge_nodes[12][2] = Hex::Hex(Node* nodes[8], unsigned value) - : Cell(MshElemType::HEXAHEDRON, value) + : Cell(value) { _nodes = nodes; _neighbors = new Element*[6]; @@ -57,7 +57,7 @@ Hex::Hex(Node* nodes[8], unsigned value) } Hex::Hex(Node* n0, Node* n1, Node* n2, Node* n3, Node* n4, Node* n5, Node* n6, Node* n7, unsigned value) - : Cell(MshElemType::HEXAHEDRON, value) + : Cell(value) { _nodes = new Node*[8]; _nodes[0] = n0; @@ -75,7 +75,7 @@ Hex::Hex(Node* n0, Node* n1, Node* n2, Node* n3, Node* n4, Node* n5, Node* n6, N } Hex::Hex(const Hex &hex) - : Cell(MshElemType::HEXAHEDRON, hex.getValue()) + : Cell(hex.getValue()) { _nodes = new Node*[8]; for (unsigned i=0; i<8; i++) @@ -114,5 +114,9 @@ const Element* Hex::getFace(unsigned i) const return NULL; } +Element* Hex::clone() const +{ + return new Hex(*this); } +} // end namespace MeshLib diff --git a/MeshLib/Elements/Hex.h b/MeshLib/Elements/Hex.h index b08aff03a36fff9fdd7f2fcb37effda80303b2a0..f38b0217a60bfc63b4c01d915425395b0472ea62 100644 --- a/MeshLib/Elements/Hex.h +++ b/MeshLib/Elements/Hex.h @@ -69,6 +69,18 @@ public: /// Get the number of nodes for this element. virtual unsigned getNNodes() const { return 8; }; + /** + * Method returns the type of the element. In this case HEXAHEDRON will be returned. + * @return MshElemType::HEXAHEDRON + */ + virtual MshElemType::type getType() const { return MshElemType::HEXAHEDRON; } + + /** + * Method clone is inherited from class Element. It makes a deep copy of the Hex instance. + * @return an exact copy of the object + */ + virtual Element* clone() const; + protected: /// Calculates the volume of a convex hexahedron by partitioning it into six tetrahedra. double computeVolume(); diff --git a/MeshLib/Elements/Prism.cpp b/MeshLib/Elements/Prism.cpp index ba84cbaedb8824062c171d9d75e8980396930f0f..ad5d9a6f8ef0ab3fc4bc0b730917ecd394588197 100644 --- a/MeshLib/Elements/Prism.cpp +++ b/MeshLib/Elements/Prism.cpp @@ -45,7 +45,7 @@ const unsigned Prism::_n_face_nodes[5] = { 3, 4, 4, 4, 3 }; Prism::Prism(Node* nodes[6], unsigned value) - : Cell(MshElemType::PRISM, value) + : Cell(value) { _nodes = nodes; _neighbors = new Element*[5]; @@ -55,7 +55,7 @@ Prism::Prism(Node* nodes[6], unsigned value) } Prism::Prism(Node* n0, Node* n1, Node* n2, Node* n3, Node* n4, Node* n5, unsigned value) - : Cell(MshElemType::PRISM, value) + : Cell(value) { _nodes = new Node*[6]; _nodes[0] = n0; @@ -71,7 +71,7 @@ Prism::Prism(Node* n0, Node* n1, Node* n2, Node* n3, Node* n4, Node* n5, unsigne } Prism::Prism(const Prism &prism) - : Cell(MshElemType::PRISM, prism.getValue()) + : Cell(prism.getValue()) { _nodes = new Node*[6]; for (unsigned i=0; i<6; i++) @@ -119,5 +119,10 @@ unsigned Prism::getNFaceNodes(unsigned i) const return 0; } +Element* Prism::clone() const +{ + return new Prism(*this); } +} // end namespace MeshLib + diff --git a/MeshLib/Elements/Prism.h b/MeshLib/Elements/Prism.h index 4397ac09532df7ce358b0aceb657e1c1410e87e8..9ae85f7c102604396571ccac2cf4ec6d01f8e6c6 100644 --- a/MeshLib/Elements/Prism.h +++ b/MeshLib/Elements/Prism.h @@ -67,6 +67,19 @@ public: /// Get the number of nodes for this element. virtual unsigned getNNodes() const { return 6; }; + /** + * Method returns the type of the element. In this case PRISM will be returned. + * @return MshElemType::PRISM + */ + virtual MshElemType::type getType() const { return MshElemType::PRISM; } + + /** + * Method clone is inherited from class Element. It makes a deep copy of the + * Hex instance employing the copy constructor of class Prism. + * @return an exact copy of the object + */ + virtual Element* clone() const; + protected: /// Calculates the volume of a prism by subdividing it into three tetrahedra. double computeVolume(); diff --git a/MeshLib/Elements/Pyramid.cpp b/MeshLib/Elements/Pyramid.cpp index def1d878c6dba8a9fa57320d5eccdab68f3eb07b..5cf90b0b1d74ca1cc09aa52d31f4ed6e452a29d7 100644 --- a/MeshLib/Elements/Pyramid.cpp +++ b/MeshLib/Elements/Pyramid.cpp @@ -44,7 +44,7 @@ const unsigned Pyramid::_n_face_nodes[5] = { 3, 3, 3, 3, 4 }; Pyramid::Pyramid(Node* nodes[5], unsigned value) - : Cell(MshElemType::PYRAMID, value) + : Cell(value) { _nodes = nodes; _neighbors = new Element*[5]; @@ -54,7 +54,7 @@ Pyramid::Pyramid(Node* nodes[5], unsigned value) } Pyramid::Pyramid(Node* n0, Node* n1, Node* n2, Node* n3, Node* n4, unsigned value) - : Cell(MshElemType::PYRAMID, value) + : Cell(value) { _nodes = new Node*[5]; _nodes[0] = n0; @@ -70,7 +70,7 @@ Pyramid::Pyramid(Node* n0, Node* n1, Node* n2, Node* n3, Node* n4, unsigned valu } Pyramid::Pyramid(const Pyramid &pyramid) - : Cell(MshElemType::PYRAMID, pyramid.getValue()) + : Cell(pyramid.getValue()) { _nodes = new Node*[5]; _neighbors = new Element*[5]; @@ -118,5 +118,9 @@ unsigned Pyramid::getNFaceNodes(unsigned i) const return 0; } +Element* Pyramid::clone() const +{ + return new Pyramid(*this); } +} // end namespace MeshLib diff --git a/MeshLib/Elements/Pyramid.h b/MeshLib/Elements/Pyramid.h index dff51812de8c2140948a3ef0c2a0f74955af0cf2..0564fff185b686a23f483fb4bce82495e607caf6 100644 --- a/MeshLib/Elements/Pyramid.h +++ b/MeshLib/Elements/Pyramid.h @@ -67,6 +67,19 @@ public: /// Get the number of nodes for this element. virtual unsigned getNNodes() const { return 5; }; + /** + * Method returns the type of the element. In this case PYRAMID will be returned. + * @return MshElemType::PYRAMID + */ + virtual MshElemType::type getType() const { return MshElemType::PYRAMID; } + + /** + * Method clone is inherited from class Element. It makes a deep copy of the + * Pyramid instance employing the copy constructor of class Pyramid. + * @return an exact copy of the object + */ + virtual Element* clone() const; + protected: /// Calculates the volume of a prism by subdividing it into two tetrahedra. double computeVolume(); diff --git a/MeshLib/Elements/Quad.cpp b/MeshLib/Elements/Quad.cpp index 385ebdb86f5faf7bc6cc170075ce73adcf733288..a4b632f16ea245b8eaa3b515c2b57f6b8592c1ce 100644 --- a/MeshLib/Elements/Quad.cpp +++ b/MeshLib/Elements/Quad.cpp @@ -28,7 +28,7 @@ const unsigned Quad::_edge_nodes[4][2] = Quad::Quad(Node* nodes[4], unsigned value) - : Face(MshElemType::TRIANGLE, value) + : Face(value) { _nodes = nodes; _neighbors = new Element*[4]; @@ -38,7 +38,7 @@ Quad::Quad(Node* nodes[4], unsigned value) } Quad::Quad(Node* n0, Node* n1, Node* n2, Node* n3, unsigned value) - : Face(MshElemType::TRIANGLE, value) + : Face(value) { _nodes = new Node*[4]; _nodes[0] = n0; @@ -52,7 +52,7 @@ Quad::Quad(Node* n0, Node* n1, Node* n2, Node* n3, unsigned value) } Quad::Quad(const Quad &quad) - : Face(MshElemType::QUAD, quad.getValue()) + : Face(quad.getValue()) { _nodes = new Node*[4]; _neighbors = new Element*[4]; @@ -74,5 +74,11 @@ double Quad::computeArea() + MathLib::calcTriangleArea(_nodes[2]->getCoords(), _nodes[3]->getCoords(), _nodes[0]->getCoords()); } +Element* Quad::clone() const +{ + return new Quad(*this); +} + + } diff --git a/MeshLib/Elements/Quad.h b/MeshLib/Elements/Quad.h index bc8069018cf619c4057dcd2398ba3542d4dd8efc..e633c5df7faf73e780db4dd31c13f043e76021dd 100644 --- a/MeshLib/Elements/Quad.h +++ b/MeshLib/Elements/Quad.h @@ -57,6 +57,18 @@ public: /// Get the number of nodes for this element. virtual unsigned getNNodes() const { return 4; }; + /** + * Method returns the type of the element. In this case QUAD will be returned. + * @return MshElemType::QUAD + */ + virtual MshElemType::type getType() const { return MshElemType::QUAD; } + + /** + * Method clone is inherited from class Element. It makes a deep copy of the Quad instance. + * @return an exact copy of the object + */ + virtual Element* clone() const; + protected: /// Calculates the area of a convex quadliteral by dividing it into two triangles. double computeArea(); diff --git a/MeshLib/Elements/Tet.cpp b/MeshLib/Elements/Tet.cpp index 13a03678d2b2ec77f07ec699b31d960c50df3d3c..d4c1aa3ef910065ce9629fd9460c0962dd167cc3 100644 --- a/MeshLib/Elements/Tet.cpp +++ b/MeshLib/Elements/Tet.cpp @@ -39,7 +39,7 @@ const unsigned Tet::_edge_nodes[6][2] = Tet::Tet(Node* nodes[4], unsigned value) - : Cell(MshElemType::TETRAHEDRON, value) + : Cell(value) { _nodes = nodes; _neighbors = new Element*[4]; @@ -49,7 +49,7 @@ Tet::Tet(Node* nodes[4], unsigned value) } Tet::Tet(Node* n0, Node* n1, Node* n2, Node* n3, unsigned value) - : Cell(MshElemType::TETRAHEDRON, value) + : Cell(value) { _nodes = new Node*[4]; _nodes[0] = n0; @@ -63,7 +63,7 @@ Tet::Tet(Node* n0, Node* n1, Node* n2, Node* n3, unsigned value) } Tet::Tet(unsigned value) - : Cell(MshElemType::TETRAHEDRON, value) + : Cell(value) { _neighbors = new Element*[4]; for (unsigned i=0; i<4; i++) @@ -71,7 +71,7 @@ Tet::Tet(unsigned value) } Tet::Tet(const Tet &tet) - : Cell(MshElemType::TETRAHEDRON, tet.getValue()) + : Cell(tet.getValue()) { _nodes = new Node*[4]; _neighbors = new Element*[4]; @@ -106,5 +106,11 @@ const Element* Tet::getFace(unsigned i) const return NULL; } +Element* Tet::clone() const +{ + return new Tet(*this); +} + + } diff --git a/MeshLib/Elements/Tet.h b/MeshLib/Elements/Tet.h index 88fd3dfbaec18a766366f3f9e93552fc26b7da77..7164b02a93a2b17e1e5a602155c4ec36b2608ea4 100644 --- a/MeshLib/Elements/Tet.h +++ b/MeshLib/Elements/Tet.h @@ -68,6 +68,18 @@ public: /// Get the number of nodes for this element. virtual unsigned getNNodes() const { return 4; }; + /** + * Method returns the type of the element. In this case TETRAHEDRON will be returned. + * @return MshElemType::TETRAHEDRON + */ + virtual MshElemType::type getType() const { return MshElemType::TETRAHEDRON; } + + /** + * Method clone is inherited from class Element. It makes a deep copy of the Tet instance. + * @return an exact copy of the object + */ + virtual Element* clone() const; + protected: /// Constructor without nodes (for use of derived classes) Tet(unsigned value = 0); diff --git a/MeshLib/Elements/Tet10.cpp b/MeshLib/Elements/Tet10.cpp index 1f7aeba54e6c40e32cd58e7225823ccfa6604fa7..55e8fa6029a173dc90a48a5449bfd6ae17da1686 100644 --- a/MeshLib/Elements/Tet10.cpp +++ b/MeshLib/Elements/Tet10.cpp @@ -59,5 +59,10 @@ void Tet10::calcCentroid() _centroid = 0; } +Element* Tet10::clone() const +{ + return new Tet10(*this); +} + } diff --git a/MeshLib/Elements/Tet10.h b/MeshLib/Elements/Tet10.h index 2f7d84ca09d4abdd62f5776753cc7eb42639fc11..f2928c9310cbb1c329ca22e34d34683570311db7 100644 --- a/MeshLib/Elements/Tet10.h +++ b/MeshLib/Elements/Tet10.h @@ -61,6 +61,18 @@ public: /// Get the number of nodes for this element. unsigned getNNodes() const { return 10; }; + /** + * Method returns the type of the element. In this case TETRAHEDRON will be returned. + * @return MshElemType::TETRAHEDRON + */ + virtual MshElemType::type getType() const { return MshElemType::TETRAHEDRON; } + + /** + * Method clone is inherited from class Element. It makes a deep copy of the Tet10 instance. + * @return an exact copy of the object + */ + virtual Element* clone() const; + protected: /// Calculates the volume of a tetrahedron via the determinant of the matrix given by its four points. void calcCentroid(); diff --git a/MeshLib/Elements/Tri.cpp b/MeshLib/Elements/Tri.cpp index 30b6e5cb4ee963b65d1782a3dbcdba33dc9b9786..3778b7c0e160609ddc2577364d04f7c443b731b0 100644 --- a/MeshLib/Elements/Tri.cpp +++ b/MeshLib/Elements/Tri.cpp @@ -27,7 +27,7 @@ const unsigned Tri::_edge_nodes[3][2] = Tri::Tri(Node* nodes[3], unsigned value) - : Face(MshElemType::TRIANGLE, value) + : Face(value) { _nodes = nodes; _neighbors = new Element*[3]; @@ -37,7 +37,7 @@ Tri::Tri(Node* nodes[3], unsigned value) } Tri::Tri(Node* n0, Node* n1, Node* n2, unsigned value) - : Face(MshElemType::TRIANGLE, value) + : Face(value) { _nodes = new Node*[3]; _nodes[0] = n0; @@ -50,7 +50,7 @@ Tri::Tri(Node* n0, Node* n1, Node* n2, unsigned value) } Tri::Tri(const Tri &tri) - : Face(MshElemType::TRIANGLE, tri.getValue()) + : Face(tri.getValue()) { _nodes = new Node*[3]; _neighbors = new Element*[3]; @@ -66,6 +66,11 @@ Tri::~Tri() { } +Element* Tri::clone() const +{ + return new Tri(*this); +} + double Tri::computeArea() { return MathLib::calcTriangleArea(_nodes[0]->getCoords(), _nodes[1]->getCoords(), _nodes[2]->getCoords()); diff --git a/MeshLib/Elements/Tri.h b/MeshLib/Elements/Tri.h index 2bb6d27397b3ddaa2eb9b1c5f953caeef320d321..3440323edb945e5c9068379babfc1a6fda810f18 100644 --- a/MeshLib/Elements/Tri.h +++ b/MeshLib/Elements/Tri.h @@ -57,6 +57,18 @@ public: /// Get the number of nodes for this element. virtual unsigned getNNodes() const { return 3; }; + /** + * Method returns the type of the element. In this case TRIANGLE will be returned. + * @return MshElemType::TRIANGLE + */ + virtual MshElemType::type getType() const { return MshElemType::TRIANGLE; } + + /** + * Method clone is inherited from class Element. It makes a deep copy of the Tri instance. + * @return an exact copy of the object + */ + virtual Element* clone() const; + protected: /// Calculates the area of the triangle by returning half of the area of the corresponding parallelogram. double computeArea(); diff --git a/MeshLib/Mesh.cpp b/MeshLib/Mesh.cpp index 98573cab59f0e7db186a6830258e4d42cb193a21..09ee63e4d3271972c6499151a35f14a637cf99f0 100644 --- a/MeshLib/Mesh.cpp +++ b/MeshLib/Mesh.cpp @@ -75,7 +75,7 @@ void Mesh::addElement(Element* elem) { _elements.push_back(elem); - // add element informatin to nodes + // add element information to nodes unsigned nNodes (elem->getNNodes()); for (unsigned i=0; i<nNodes; i++) elem->_nodes[i]->addElement(elem); diff --git a/MeshLib/MeshCoarsener.cpp b/MeshLib/MeshCoarsener.cpp new file mode 100644 index 0000000000000000000000000000000000000000..8abd103d3f7b1320c5a9649247b8a578d02460e0 --- /dev/null +++ b/MeshLib/MeshCoarsener.cpp @@ -0,0 +1,127 @@ +/** + * Copyright (c) 2012, OpenGeoSys Community (http://www.opengeosys.net) + * Distributed under a Modified BSD License. + * See accompanying file LICENSE.txt or + * http://www.opengeosys.net/LICENSE.txt + * + * \file MeshCoarsener.cpp + * + * Created on Aug 3, 2012 by Thomas Fischer + */ + +#include "MeshCoarsener.h" + +// GeoLib +#include "Grid.h" + +// MeshLib +#include "Mesh.h" +#include "Node.h" +#include "Elements/Element.h" + +namespace MeshLib { + +MeshCoarsener::MeshCoarsener(Mesh const*const mesh) : + _orig_mesh (mesh) +{} + +MeshCoarsener::~MeshCoarsener() +{} + +Mesh* MeshCoarsener::operator()(double min_distance) +{ + // copy mesh nodes, reset mesh node ids and + // store original mesh node ids in a vector + std::vector<Node*> const& orig_nodes(_orig_mesh->getNodes()); + const size_t n_nodes(orig_nodes.size()); + std::vector<Node*> nodes(n_nodes); + std::map<size_t,size_t> orig_ids_map; + for (size_t k(0); k < n_nodes; k++) { + nodes[k] = new Node(orig_nodes[k]->getCoords(), k); + orig_ids_map.insert(std::pair<size_t, size_t>(orig_nodes[k]->getID(), k)); + } + + // init grid + GeoLib::Grid<Node>* grid(new GeoLib::Grid<Node>(nodes, 64)); + + // init id map + std::vector<size_t> id_map(n_nodes); + for (size_t k(0); k < n_nodes; k++) { + id_map[k] = nodes[k]->getID(); + } + + const double sqr_min_distance (min_distance * min_distance); + + // do the work - search nearest nodes + for (size_t k(0); k < n_nodes; k++) { + std::vector<std::vector<Node*> const*> node_vecs_intersecting_cube; + Node const*const node(nodes[k]); + const size_t node_id(node->getID()); + grid->getVecsOfGridCellsIntersectingCube(node->getCoords(), min_distance, node_vecs_intersecting_cube); + + const size_t n_vecs (node_vecs_intersecting_cube.size()); + for (size_t i(0); i<n_vecs; i++) { + std::vector<Node*> const* node_vec (node_vecs_intersecting_cube[i]); + const size_t n_loc_nodes (node_vec->size()); + for (size_t j(0); j<n_loc_nodes; j++) { + if (node_id < (*node_vec)[j]->getID()) { + if (MathLib::sqrDist(node->getCoords(), (*node_vec)[j]->getCoords()) < sqr_min_distance) { + // two nodes are very close to each other + id_map[j] = k; + } + } + } + } + } + + // apply changes to id_map + for (size_t k(0), cnt(0); k < n_nodes; k++) { + if (id_map[k] != k) { + if (id_map[k] != id_map[id_map[k]]) { + id_map[k] = id_map[id_map[k]]; + } + } else { + id_map[k] = cnt++; + } + } + + // delete unused nodes + for (size_t k(0); k < n_nodes; k++) { + if (id_map[k] != k) { + delete nodes[k]; + nodes[k] = NULL; + } + } + + // remove NULL-ptr from node vector + for (std::vector<Node*>::iterator it(nodes.begin()); it != nodes.end(); ) { + if (*it == NULL) { + it = nodes.erase (it); + } else { + it++; + } + } + + // copy mesh elements, reset the node pointers + std::vector<Element*> const& orig_elements(_orig_mesh->getElements()); + const size_t n_elements(orig_elements.size()); + std::vector<Element*> elements(n_elements); + for (size_t k(0); k < n_elements; k++) { + elements[k] = orig_elements[k]->clone(); + const size_t n_nodes_element (elements[k]->getNNodes()); + for (size_t i(0); i<n_nodes_element; i++) { + const size_t orig_node_id (elements[k]->getNode(i)->getID()); + size_t orig_node_pos (std::numeric_limits<size_t>::max()); + std::map<size_t, size_t>::const_iterator it(orig_ids_map.find(orig_node_id)); + if (it == orig_ids_map.end()) { + std::cerr << "[MeshCoarsener::operator()] could not found mesh node id" << std::endl; + } + orig_node_pos = it->second; + elements[k]->setNode(i, nodes[id_map[orig_node_pos]]); + } + } + + return new Mesh ("test", nodes, elements); +} + +} // end namespace MeshLib diff --git a/MeshLib/MeshCoarsener.h b/MeshLib/MeshCoarsener.h new file mode 100644 index 0000000000000000000000000000000000000000..4704f3c74946e593a3f525f2070f7610f1012e4a --- /dev/null +++ b/MeshLib/MeshCoarsener.h @@ -0,0 +1,54 @@ +/** + * Copyright (c) 2012, OpenGeoSys Community (http://www.opengeosys.net) + * Distributed under a Modified BSD License. + * See accompanying file LICENSE.txt or + * http://www.opengeosys.net/LICENSE.txt + * + * \file MeshCoarsener.h + * + * Created on Aug 3, 2012 by Thomas Fischer + */ + +#ifndef MESHCOARSENER_H_ +#define MESHCOARSENER_H_ + +#include <vector> +#include <cstddef> + +// forward declaration +namespace MeshLib { +class Mesh; +} + +namespace MeshLib { + +/** + * The class MeshCoarsener merges mesh nodes that have a smaller + * distance than a (user) given minimal distance. + */ +class MeshCoarsener { +public: + /** + * Constructor of class MeshCoarsener that takes the mesh object that should be coarsened. + * @param mesh the mesh object + */ + MeshCoarsener(Mesh const*const mesh); + + /** + * destructor. + */ + virtual ~MeshCoarsener(); + + /** + * create new mesh and apply the coarsening process to the mesh + * @param min_distance + */ + Mesh* operator() (double min_distance); + +private: + Mesh const*const _orig_mesh; +}; + +} + +#endif /* MESHCOARSENER_H_ */ diff --git a/SimpleTests/MeshTests/CMakeLists.txt b/SimpleTests/MeshTests/CMakeLists.txt index 1da4fad39113a8053f9615ca37155fd2c7ed9b30..1dcf599294e01ce22b2da7aee589d45cc32bac6f 100644 --- a/SimpleTests/MeshTests/CMakeLists.txt +++ b/SimpleTests/MeshTests/CMakeLists.txt @@ -25,3 +25,19 @@ TARGET_LINK_LIBRARIES ( MeshRead logog ) +# Create CollapseMeshNodes executable +ADD_EXECUTABLE( CollapseMeshNodes + CollapseMeshNodes.cpp + ${SOURCES} + ${HEADERS} +) + +TARGET_LINK_LIBRARIES ( CollapseMeshNodes + MeshLib + FileIO + MathLib + BaseLib + GeoLib + logog +) + diff --git a/SimpleTests/MeshTests/CollapseMeshNodes.cpp b/SimpleTests/MeshTests/CollapseMeshNodes.cpp new file mode 100644 index 0000000000000000000000000000000000000000..add35350f6107dc262927a32c074434913c4cf1a --- /dev/null +++ b/SimpleTests/MeshTests/CollapseMeshNodes.cpp @@ -0,0 +1,92 @@ +/** + * Copyright (c) 2012, OpenGeoSys Community (http://www.opengeosys.net) + * Distributed under a Modified BSD License. + * See accompanying file LICENSE.txt or + * http://www.opengeosys.net/LICENSE.txt + * + * \file CollapseMeshNodes.cpp + * + * Created on Jul 31, 2012 by Thomas Fischer + */ + +// BaseLib +#include "MemWatch.h" +#include "RunTime.h" +#include "tclap/CmdLine.h" + +// BaseLib/logog +#include "logog.hpp" + + +// MeshLib +#include "Node.h" +#include "Elements/Element.h" +#include "Mesh.h" +#include "MeshIO.h" +#include "MeshCoarsener.h" + +int main(int argc, char *argv[]) +{ + LOGOG_INITIALIZE(); + logog::Cout* logogCout = new logog::Cout; + + TCLAP::CmdLine cmd("Collapse mesh nodes and, if necessary, remove elements", ' ', "0.1"); + + // Define a value argument and add it to the command line. + // A value arg defines a flag and a type of value that it expects, + // such as "-m meshfile". + TCLAP::ValueArg<std::string> input_mesh_arg("m","mesh","input mesh file",true,"","string"); + + // Add the argument mesh_arg to the CmdLine object. The CmdLine object + // uses this Arg to parse the command line. + cmd.add( input_mesh_arg ); + + cmd.parse( argc, argv ); + + std::string fname (input_mesh_arg.getValue()); + + FileIO::MeshIO mesh_io; +#ifndef WIN32 + BaseLib::MemWatch mem_watch; + unsigned long mem_without_mesh (mem_watch.getVirtMemUsage()); + BaseLib::RunTime run_time; + run_time.start(); +#endif + MeshLib::Mesh* mesh = mesh_io.loadMeshFromFile(fname); +#ifndef WIN32 + if (mesh) { + unsigned long mem_with_mesh (mem_watch.getVirtMemUsage()); +// std::cout << "mem for mesh: " << (mem_with_mesh - mem_without_mesh)/(1024*1024) << " MB" << std::endl; + INFO ("mem for mesh: %i MB", (mem_with_mesh - mem_without_mesh)/(1024*1024)); + } + run_time.stop(); +// std::cout << "time for reading: " << run_time.elapsed() << " s" << std::endl; + if (mesh) { + INFO ("time for reading: %f s", run_time.elapsed()); + } +#endif + +#ifndef WIN32 + unsigned long mem_without_meshgrid (mem_watch.getVirtMemUsage()); + run_time.start(); +#endif + MeshLib::MeshCoarsener mesh_coarsener(mesh); + MeshLib::Mesh *collapsed_mesh(mesh_coarsener (10)); +#ifndef WIN32 + run_time.stop(); + unsigned long mem_with_meshgrid (mem_watch.getVirtMemUsage()); + INFO ("mem for meshgrid: %i MB", (mem_with_meshgrid - mem_without_meshgrid)/(1024*1024)); + INFO ("time for collapsing: %f s", run_time.elapsed()); +#endif + + mesh_io.setMesh(collapsed_mesh); + std::string out_fname("/home/fischeth/workspace/OGS-6/Build/CollapsedMesh.msh"); + INFO ("writing collapsed mesh to %s", out_fname.c_str()); + mesh_io.writeToFile(out_fname); + INFO ("done"); + + delete mesh; + delete collapsed_mesh; + delete logogCout; + LOGOG_SHUTDOWN(); +}