diff --git a/FileIO/Legacy/MeshIO.cpp b/FileIO/Legacy/MeshIO.cpp
index 3c047059cfee40f08c88b4c90f05e6b4039feca3..5aa53bd677eb9761c9e3f3d3c4593fc96242ad84 100644
--- a/FileIO/Legacy/MeshIO.cpp
+++ b/FileIO/Legacy/MeshIO.cpp
@@ -231,14 +231,11 @@ void MeshIO::writeElementsExceptLines(std::vector<MeshLib::Element*> const& ele_
}
}
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()) << " ";
- const size_t nElemNodes (ele_vec[i]->getNNodes()-1);
- for(size_t j = 0; j < nElemNodes; j++)
- out << ele_vec[i]->getNode(nElemNodes-j-1)->getID() << " ";
+ for (size_t i(0), k(0); i < ele_vector_size; i++) {
+ if (non_line_element[i] && non_null_element[i]) {
+ out << k << " " << ele_vec[i]->getValue() << " " << 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++;
}
diff --git a/MathLib/MathTools.h b/MathLib/MathTools.h
index 505e58f990014a9abe05c417b2c5ad71900bd168..ad83ef5618b9aecbe10617b0cf21176a21d33154 100644
--- a/MathLib/MathTools.h
+++ b/MathLib/MathTools.h
@@ -121,7 +121,7 @@ double sqrDist(const double* p0, const double* p1);
float normalize(float min, float max, float val);
/**
- * Let \f$p_0, p_1, p_2 \in \mathbb{R}^3\f$. The function getAngle
+ * Let \f$p_0, p_1, p_2 \in R^3\f$. The function getAngle
* computes the angle between the edges \f$(p_0,p_1)\f$ and \f$(p_1,p_2)\f$
* @param p0 start point of edge 0
* @param p1 end point of edge 0 and start point of edge 1
diff --git a/MeshLib/CMakeLists.txt b/MeshLib/CMakeLists.txt
index 537a0e1d6ad731059adac494570d900211dc4660..51087b6a8c68ab5ddbd7f07ebf06ff089938f8d0 100644
--- a/MeshLib/CMakeLists.txt
+++ b/MeshLib/CMakeLists.txt
@@ -12,6 +12,7 @@ include_directories(
../BaseLib
../GeoLib
../MathLib
+ ${CMAKE_SOURCE_DIR}/BaseLib/logog/include
)
@@ -19,5 +20,6 @@ target_link_libraries (MeshLib
BaseLib
GeoLib
MathLib
+ logog
)
diff --git a/MeshLib/Elements/Edge.cpp b/MeshLib/Elements/Edge.cpp
index 3c0dca2874215119989f50817311895a42a93a96..a5db2266f985feb9e6bb1d0d53a4488d481f4139 100644
--- a/MeshLib/Elements/Edge.cpp
+++ b/MeshLib/Elements/Edge.cpp
@@ -57,5 +57,14 @@ Element* Edge::clone() const
return new Edge(*this);
}
+Element* Edge::reviseElement() const
+{
+ if (_nodes[0] == _nodes[1]) {
+ return NULL;
+ }
+
+ return NULL;
+}
+
}
diff --git a/MeshLib/Elements/Edge.h b/MeshLib/Elements/Edge.h
index cc6894d32cd0341fabd3849dcadeba95f7420963..5a4461d0b7a7576bf305fc223412f91f2e537fd0 100644
--- a/MeshLib/Elements/Edge.h
+++ b/MeshLib/Elements/Edge.h
@@ -22,10 +22,7 @@ class Node;
/**
* A 1d Edge or Line Element.
* @code
- *
- * Edge: o--------o
- * 0 1
- *
+ * 0--------1
* @endcode
*/
class Edge : public Element
@@ -80,12 +77,18 @@ public:
virtual Element* clone() const;
+ /**
+ * If for instance two nodes of the element are collapsed the Edge element disappears.
+ * @return NULL
+ */
+ virtual Element* reviseElement() const;
+
protected:
/// Calculate the length of this 1d element.
double computeLength();
/// 1D elements have no edges.
- Node* getEdgeNode(unsigned edge_id, unsigned node_id) const { (void)edge_id; (void)node_id; return NULL; };
+ Node const* getEdgeNode(unsigned edge_id, unsigned node_id) const { (void)edge_id; (void)node_id; return NULL; };
/// 1D elements have no faces.
Node* getFaceNode(unsigned face_id, unsigned node_id) const { (void)face_id; (void)node_id; return NULL; };
diff --git a/MeshLib/Elements/Element.cpp b/MeshLib/Elements/Element.cpp
index a8572c56414295b49e0cd7917363f116ad5a3588..c751a180e46254e6c197d678022153be1c142ee9 100644
--- a/MeshLib/Elements/Element.cpp
+++ b/MeshLib/Elements/Element.cpp
@@ -72,8 +72,8 @@ const Element* Element::getEdge(unsigned i) const
if (i < getNEdges())
{
Node** nodes = new Node*[2];
- nodes[0] = getEdgeNode(i,0);
- nodes[1] = getEdgeNode(i,1);
+ nodes[0] = const_cast<Node*>(getEdgeNode(i,0));
+ nodes[1] = const_cast<Node*>(getEdgeNode(i,1));
return new Edge(nodes);
}
std::cerr << "Error in MeshLib::Element::getEdge() - Index does not exist." << std::endl;
@@ -87,7 +87,7 @@ void Element::computeSqrEdgeLengthRange(double &min, double &max) const
unsigned nEdges (this->getNEdges());
for (unsigned i=0; i<nEdges; i++)
{
- double dist (MathLib::sqrDist(getEdgeNode(i,0)->getCoords(), getEdgeNode(i,1)->getCoords()));
+ double dist (MathLib::sqrDist(getEdgeNode(i,0), getEdgeNode(i,1)));
min = (dist<min) ? dist : min;
max = (dist>max) ? dist : max;
}
diff --git a/MeshLib/Elements/Element.h b/MeshLib/Elements/Element.h
index 3913348f83dff5ee063d79f45b104d90470aa7ef..20dd2b8b994fac02e3a9e4a95ba5e09eaa6b3481 100644
--- a/MeshLib/Elements/Element.h
+++ b/MeshLib/Elements/Element.h
@@ -36,16 +36,24 @@ public:
virtual void computeSqrEdgeLengthRange(double &min, double &max) const;
/**
- * \brief Tries to add an element e as neighbour to this element.
+ * \brief Tries to add an element e as neighbour to this element.
* If the elements really are neighbours, the element is added to the
- * neighbour-ist and true is returned. Otherwise false is returned.
+ * neighbour-list and true is returned. Otherwise false is returned.
*/
bool addNeighbor(Element* e);
/// Returns the length, area or volume of a 1D, 2D or 3D element
virtual double getContent() const = 0;
- /// Get node with local index i.
+ /**
+ * Get node with local index i where i should be at most the number
+ * of nodes of the element
+ * @param i local index of node, at most the number of nodes of the
+ * element that you can obtain with Element::getNNodes()
+ * @return a pointer to the appropriate (and constant, i.e. not
+ * modifiable by the user) instance of class Node or a NULL pointer
+ * @sa Element::getNodeIndex()
+ */
const Node* getNode(unsigned i) const;
/**
@@ -85,7 +93,14 @@ public:
/// Get the number of nodes for this element.
virtual unsigned getNNodes() const = 0;
- /// Get the global index for the node with local index i.
+ /**
+ * Get the global index for the Node with local index i.
+ * The index i should be at most the number of nodes of the element.
+ * @param i local index of Node, at most the number of nodes of the
+ * element that you can obtain with Element::getNNodes()
+ * @return the global index or std::numeric_limits<unsigned>::max()
+ * @sa Element::getNode()
+ */
unsigned getNodeIndex(unsigned i) const;
/// Get the type of the mesh element (as a MshElemType-enum).
@@ -106,12 +121,22 @@ public:
*/
virtual Element* clone() const = 0;
+ /**
+ * This method should be called after at least two nodes of an element
+ * are collapsed. The node collapsing can/have to lead to an edge collapse.
+ * This method tries to create a new element of an appropriate type. The
+ * value of the attribute _value is carried over. In contrast to this the
+ * neighbor information is not carried over.
+ * @return an element of a different element type (MshElemType) or NULL
+ */
+ virtual Element* reviseElement() const = 0;
+
protected:
/// Constructor for a generic mesh element without an array of mesh nodes.
Element(unsigned value = 0);
/// Return a specific edge node.
- virtual Node* getEdgeNode(unsigned edge_id, unsigned node_id) const = 0;
+ virtual Node const* getEdgeNode(unsigned edge_id, unsigned node_id) const = 0;
Node** _nodes;
unsigned _value;
diff --git a/MeshLib/Elements/Hex.cpp b/MeshLib/Elements/Hex.cpp
index f296292fd41373c18448b7070fee3a2cf497b696..b12fd76501bb10cf3c3e5f204997ad8b2f2da2de 100644
--- a/MeshLib/Elements/Hex.cpp
+++ b/MeshLib/Elements/Hex.cpp
@@ -13,6 +13,7 @@
#include "Hex.h"
#include "Node.h"
#include "Quad.h"
+#include "Prism.h"
#include "MathTools.h"
@@ -119,4 +120,50 @@ Element* Hex::clone() const
return new Hex(*this);
}
+Element* Hex::reviseElement() const
+{
+ std::vector<size_t> collapsed_edges;
+ for (size_t edge(0); edge<getNEdges(); edge++) {
+ if (_nodes[_edge_nodes[edge][0]] == _nodes[_edge_nodes[edge][1]]) {
+ collapsed_edges.push_back(edge);
+ }
+ }
+
+ if (collapsed_edges.size() == 1) {
+ std::cerr << "[Hex::reviseElement()] collapsing of one edge in hexahedron not handled" << std::endl;
+ return NULL;
+ }
+
+ if (collapsed_edges.size() == 2) {
+ // try to create a prism out of the hex
+ if (collapsed_edges[0] == 0 && collapsed_edges[1] == 2) {
+ return new Prism(_nodes[0],_nodes[4], _nodes[5], _nodes[3], _nodes[7], _nodes[6], _value);
+ }
+ if (collapsed_edges[0] == 1 && collapsed_edges[1] == 3) {
+ return new Prism(_nodes[0],_nodes[4], _nodes[7], _nodes[1], _nodes[5], _nodes[6], _value);
+ }
+ if (collapsed_edges[0] == 4 && collapsed_edges[1] == 5) {
+ return new Prism(_nodes[0],_nodes[7], _nodes[3], _nodes[1], _nodes[6], _nodes[2], _value);
+ }
+ if (collapsed_edges[0] == 5 && collapsed_edges[1] == 6) {
+ return new Prism(_nodes[0],_nodes[1], _nodes[4], _nodes[3], _nodes[2], _nodes[7], _value);
+ }
+ if (collapsed_edges[0] == 6 && collapsed_edges[1] == 7) {
+ return new Prism(_nodes[0],_nodes[3], _nodes[4], _nodes[1], _nodes[2], _nodes[5], _value);
+ }
+ if (collapsed_edges[0] == 7 && collapsed_edges[1] == 4) {
+ return new Prism(_nodes[0],_nodes[1], _nodes[5], _nodes[3], _nodes[2], _nodes[6], _value);
+ }
+ if (collapsed_edges[0] == 8 && collapsed_edges[1] == 10) {
+ return new Prism(_nodes[0],_nodes[1], _nodes[4], _nodes[3], _nodes[2], _nodes[7], _value);
+ }
+ if (collapsed_edges[0] == 9 && collapsed_edges[1] == 11) {
+ return new Prism(_nodes[0],_nodes[3], _nodes[4], _nodes[1], _nodes[2], _nodes[5], _value);
+ }
+ return NULL;
+ }
+
+ return NULL;
+}
+
} // end namespace MeshLib
diff --git a/MeshLib/Elements/Hex.h b/MeshLib/Elements/Hex.h
index f38b0217a60bfc63b4c01d915425395b0472ea62..172c9200604f1ccf1f906563a687cd1415074228 100644
--- a/MeshLib/Elements/Hex.h
+++ b/MeshLib/Elements/Hex.h
@@ -21,18 +21,24 @@ namespace MeshLib {
* A 3d Hexahedron Element.
* @code
*
- * Hex: 7 6
- * o--------o
- * /: /|
- * / : / |
- * 4 / : 5 / |
- * o--------o |
- * | o....|...o 2
- * | .3 | /
- * | . | /
- * |. |/
- * o--------o
- * 0 1
+ * Hex:
+ * 10
+ * 7-----------6
+ * /: /|
+ * / : / |
+ * 11/ : /9 |
+ * / 7: / | 6
+ * / : 8 / |
+ * 4-----------5 |
+ * | : | 2 |
+ * | 3.....|.....2
+ * | . | /
+ * 4 | . |5 /
+ * | 3. | / 1
+ * | . | /
+ * |. |/
+ * 0-----------1
+ * 0
*
* @endcode
*/
@@ -81,12 +87,18 @@ public:
*/
virtual Element* clone() const;
+ /**
+ * Change the element type from hexahedron to a prism if two appropriate edges of the hexahedron are collapsed.
+ * @return a prism element with nice properties or NULL
+ */
+ virtual Element* reviseElement() const;
+
protected:
/// Calculates the volume of a convex hexahedron by partitioning it into six tetrahedra.
double computeVolume();
/// Return a specific edge node.
- inline Node* getEdgeNode(unsigned edge_id, unsigned node_id) const { return _nodes[_edge_nodes[edge_id][node_id]]; };
+ inline Node const* getEdgeNode(unsigned edge_id, unsigned node_id) const { return _nodes[_edge_nodes[edge_id][node_id]]; };
static const unsigned _face_nodes[6][4];
static const unsigned _edge_nodes[12][2];
diff --git a/MeshLib/Elements/Prism.cpp b/MeshLib/Elements/Prism.cpp
index ad5d9a6f8ef0ab3fc4bc0b730917ecd394588197..dde1cd3193d573678aa27273a59c8ad1cb91faf5 100644
--- a/MeshLib/Elements/Prism.cpp
+++ b/MeshLib/Elements/Prism.cpp
@@ -13,6 +13,7 @@
#include "Prism.h"
#include "Node.h"
#include "Tri.h"
+#include "Pyramid.h"
#include "Quad.h"
#include "MathTools.h"
@@ -124,5 +125,23 @@ Element* Prism::clone() const
return new Prism(*this);
}
+Element* Prism::reviseElement() const
+{
+ // try to create Pyramid
+ if (_nodes[_edge_nodes[3][0]] == _nodes[_edge_nodes[3][1]]) {
+ return new Pyramid(_nodes[1], _nodes[4], _nodes[5], _nodes[2], _nodes[0], _value);
+ }
+
+ if (_nodes[_edge_nodes[4][0]] == _nodes[_edge_nodes[4][1]]) {
+ return new Pyramid(_nodes[0], _nodes[2], _nodes[5], _nodes[3], _nodes[1], _value);
+ }
+
+ if (_nodes[_edge_nodes[5][0]] == _nodes[_edge_nodes[5][1]]) {
+ return new Pyramid(_nodes[0], _nodes[1], _nodes[4], _nodes[3], _nodes[2], _value);
+ }
+
+ return NULL;
+}
+
} // end namespace MeshLib
diff --git a/MeshLib/Elements/Prism.h b/MeshLib/Elements/Prism.h
index 9ae85f7c102604396571ccac2cf4ec6d01f8e6c6..7649559424aabd8c46268d1cea297938f9b13bfb 100644
--- a/MeshLib/Elements/Prism.h
+++ b/MeshLib/Elements/Prism.h
@@ -18,19 +18,25 @@
namespace MeshLib {
/**
- * A 3d Prism Element.
+ * This class represents a 3d prism element. The following sketch shows the node and edge numbering.
+ * @anchor PrismNodeAndEdgeNumbering
* @code
- *
- * Prism: 5
- * o
- * /:\
- * / : \
- * / o \
- * 3 o-------o 4
- * | . 2 . |
- * |. .|
- * o-------o
- * 0 1
+ * 5
+ * / \
+ * / : \
+ * 8/ : \7
+ * / :5 \
+ * / : 6 \
+ * 3-----------4
+ * | : |
+ * | 2 |
+ * | . . |
+ * 3| . . |4
+ * | 2. .1 |
+ * | . . |
+ * |. .|
+ * 0-----------1
+ * 0
*
* @endcode
*/
@@ -80,12 +86,23 @@ public:
*/
virtual Element* clone() const;
+ /**
+ * This method should be called after at least two nodes of the prism
+ * element are collapsed. As a consequence of the node collapsing an edge
+ * of the prism will be collapsed. If one of the edges 3, 4 or 5 (see
+ * sketch @ref PrismNodeAndEdgeNumbering) is collapsed we obtain a
+ * pyramid. In this case the method will create the appropriate
+ * object of class Pyramid.
+ * @return a pyramid object or NULL
+ */
+ virtual Element* reviseElement() const;
+
protected:
/// Calculates the volume of a prism by subdividing it into three tetrahedra.
double computeVolume();
/// Return a specific edge node.
- inline Node* getEdgeNode(unsigned edge_id, unsigned node_id) const { return _nodes[_edge_nodes[edge_id][node_id]]; };
+ inline Node const* getEdgeNode(unsigned edge_id, unsigned node_id) const { return _nodes[_edge_nodes[edge_id][node_id]]; };
static const unsigned _face_nodes[5][4];
static const unsigned _edge_nodes[9][2];
diff --git a/MeshLib/Elements/Pyramid.cpp b/MeshLib/Elements/Pyramid.cpp
index 5cf90b0b1d74ca1cc09aa52d31f4ed6e452a29d7..55b1cf0edbf07e0ff9c2e80b176c4d7cc61c2e40 100644
--- a/MeshLib/Elements/Pyramid.cpp
+++ b/MeshLib/Elements/Pyramid.cpp
@@ -13,6 +13,7 @@
#include "Pyramid.h"
#include "Node.h"
#include "Tri.h"
+#include "Tet.h"
#include "Quad.h"
#include "MathTools.h"
@@ -123,4 +124,20 @@ Element* Pyramid::clone() const
return new Pyramid(*this);
}
+Element* Pyramid::reviseElement() const
+{
+ // try to create tetrahedron
+ if (_nodes[_edge_nodes[0][0]] == _nodes[_edge_nodes[0][1]]
+ || _nodes[_edge_nodes[1][0]] == _nodes[_edge_nodes[1][1]]) {
+ return new Tet(_nodes[0], _nodes[2], _nodes[3], _nodes[4], _value);
+ }
+
+ if (_nodes[_edge_nodes[2][0]] == _nodes[_edge_nodes[2][1]]
+ || _nodes[_edge_nodes[3][0]] == _nodes[_edge_nodes[3][1]]) {
+ return new Tet(_nodes[0], _nodes[1], _nodes[2], _nodes[4], _value);
+ }
+
+ return NULL;
+}
+
} // end namespace MeshLib
diff --git a/MeshLib/Elements/Pyramid.h b/MeshLib/Elements/Pyramid.h
index 0564fff185b686a23f483fb4bce82495e607caf6..4345d8bd3d0f9c9f0e6f16aec6bad1a052001012 100644
--- a/MeshLib/Elements/Pyramid.h
+++ b/MeshLib/Elements/Pyramid.h
@@ -18,20 +18,24 @@
namespace MeshLib {
/**
- * A 3d Pyramid Element.
+ * This class represents a 3d pyramid element. The following sketch shows the node and edge numbering.
+ * @anchor PyramidNodeAndEdgeNumbering
* @code
*
- * Pyramid: o 4
- * //|\
- * // | \
- * // | \
- * 3 o/...|...o 2
- * ./ | /
- * ./ | /
- * ./ |/
- * o--------o
- * 0 1
- *
+ * 4
+ * //|\
+ * // | \
+ * 7// | \6
+ * // |5 \
+ * // | \
+ * 3/.... |.....2
+ * ./ | 2 /
+ * ./4 | /
+ * 3./ | /1
+ * ./ | /
+ * ./ |/
+ * 0------------1
+ * 0
* @endcode
*/
class Pyramid : public Cell
@@ -80,12 +84,23 @@ public:
*/
virtual Element* clone() const;
+ /**
+ * This method should be called after at least two nodes of the pyramid
+ * element are collapsed. As a consequence of the node collapsing an edge
+ * of the pyramid will be collapsed. If one of the edges 0, 1, 2 or 3 (see
+ * sketch @ref PyramidNodeAndEdgeNumbering) is collapsed we obtain a
+ * tetrahedron. In this case the method will create the appropriate
+ * object of class Tetrahedron.
+ * @return a Tetrahedron object or NULL
+ */
+ virtual Element* reviseElement() const;
+
protected:
/// Calculates the volume of a prism by subdividing it into two tetrahedra.
double computeVolume();
/// Return a specific edge node.
- inline Node* getEdgeNode(unsigned edge_id, unsigned node_id) const { return _nodes[_edge_nodes[edge_id][node_id]]; };
+ inline Node const* getEdgeNode(unsigned edge_id, unsigned node_id) const { return _nodes[_edge_nodes[edge_id][node_id]]; };
static const unsigned _face_nodes[5][4];
static const unsigned _edge_nodes[8][2];
diff --git a/MeshLib/Elements/Quad.cpp b/MeshLib/Elements/Quad.cpp
index a4b632f16ea245b8eaa3b515c2b57f6b8592c1ce..922871ddbee54da7015e137f04d08dbbc75e3991 100644
--- a/MeshLib/Elements/Quad.cpp
+++ b/MeshLib/Elements/Quad.cpp
@@ -12,6 +12,7 @@
#include "Quad.h"
#include "Node.h"
+#include "Tri.h"
#include "MathTools.h"
@@ -22,7 +23,7 @@ const unsigned Quad::_edge_nodes[4][2] =
{
{0, 1}, // Edge 0
{1, 2}, // Edge 1
- {0, 2}, // Edge 2
+ {2, 3}, // Edge 2
{0, 3} // Edge 3
};
@@ -79,6 +80,19 @@ Element* Quad::clone() const
return new Quad(*this);
}
+Element* Quad::reviseElement() const
+{
+ if (_nodes[0] == _nodes[1] || _nodes[1] == _nodes[2]) {
+ return new Tri(_nodes[0], _nodes[2], _nodes[3], _value);
+ }
+
+ if (_nodes[2] == _nodes[3] || _nodes[3] == _nodes[0]) {
+ return new Tri(_nodes[0], _nodes[1], _nodes[2], _value);
+ }
+
+ // this should not happen
+ return NULL;
+}
}
diff --git a/MeshLib/Elements/Quad.h b/MeshLib/Elements/Quad.h
index e633c5df7faf73e780db4dd31c13f043e76021dd..2dfef14ab9ad8e27589d2a809bf163fd2dce1a49 100644
--- a/MeshLib/Elements/Quad.h
+++ b/MeshLib/Elements/Quad.h
@@ -18,19 +18,18 @@
namespace MeshLib {
/**
- * A 2d Quadliteral Element.
+ * This class represents a 2d quadliteral element. The following sketch shows the node and edge numbering.
+ * @anchor QuadNodeAndEdgeNumbering
* @code
- *
- * 3 2
- * Quad: o-----------o
- * | |
+ * 2
+ * 3-----------2
* | |
* | |
+ * 3| |1
* | |
* | |
- * o-----------o
- * 0 1
- *
+ * 0-----------1
+ * 0
* @endcode
*/
class Quad : public Face
@@ -69,12 +68,23 @@ public:
*/
virtual Element* clone() const;
+ /**
+ * This method should be called after at least two nodes of the quad
+ * element are collapsed. As a consequence of the node collapsing an edge
+ * of the quad will be collapsed. If one of the edges (see
+ * sketch @ref PyramidNodeAndEdgeNumbering) is collapsed we obtain a
+ * triangle. In this case the method will create the appropriate
+ * object of class Tri.
+ * @return a Tri object or NULL
+ */
+ virtual Element* reviseElement() const;
+
protected:
/// Calculates the area of a convex quadliteral by dividing it into two triangles.
double computeArea();
/// Return a specific edge node.
- inline Node* getEdgeNode(unsigned edge_id, unsigned node_id) const { return _nodes[_edge_nodes[edge_id][node_id]]; };
+ inline Node const* getEdgeNode(unsigned edge_id, unsigned node_id) const { return _nodes[_edge_nodes[edge_id][node_id]]; };
static const unsigned _edge_nodes[4][2];
diff --git a/MeshLib/Elements/Tet.cpp b/MeshLib/Elements/Tet.cpp
index d4c1aa3ef910065ce9629fd9460c0962dd167cc3..966ec3cfb5c7522c5d7d9c11b49e8f9135aa3406 100644
--- a/MeshLib/Elements/Tet.cpp
+++ b/MeshLib/Elements/Tet.cpp
@@ -111,6 +111,23 @@ Element* Tet::clone() const
return new Tet(*this);
}
+Element* Tet::reviseElement() const
+{
+ if (_nodes[0] == _nodes[1] || _nodes[1] == _nodes[2]) {
+ return new Tri(_nodes[0], _nodes[2], _nodes[3], _value);
+ }
+
+ if (_nodes[2] == _nodes[0]) {
+ return new Tri(_nodes[0], _nodes[1], _nodes[3], _value);
+ }
+ if (_nodes[0] == _nodes[3] || _nodes[1] == _nodes[3] || _nodes[2] == _nodes[3]) {
+ return new Tri(_nodes[0], _nodes[1], _nodes[2], _value);
+ }
+
+ // this should not happen
+ return NULL;
}
+} // end namespace MeshLib
+
diff --git a/MeshLib/Elements/Tet.h b/MeshLib/Elements/Tet.h
index 7164b02a93a2b17e1e5a602155c4ec36b2608ea4..6a8ac29f00882916f8889df97ca0ac967761c31d 100644
--- a/MeshLib/Elements/Tet.h
+++ b/MeshLib/Elements/Tet.h
@@ -18,20 +18,22 @@
namespace MeshLib {
/**
- * A 3d Tetrahedron Element.
+ * This class represents a 3d tetrahedron element. The following sketch shows the node and edge numbering.
+ * @anchor TetrahedronNodeAndEdgeNumbering
* @code
- *
- * Tet: 3
- * o
- * /|\
- * / | \
- * / | \
- * 0 o...|...o 2
- * \ | /
- * \ | /
- * \|/
- * o
- * 1
+ * 3
+ * /|\
+ * / | \
+ * 3/ | \5
+ * / |4 \
+ * / | \
+ * 0.....|.....2
+ * \ | 2 /
+ * \ | /
+ * 0\ | /1
+ * \ | /
+ * \|/
+ * 1
*
* @endcode
*/
@@ -80,6 +82,15 @@ public:
*/
virtual Element* clone() const;
+ /**
+ * This method should be called after at least two nodes of the tetrahedron
+ * element are collapsed. As a consequence of the node collapsing an edge
+ * of the tetrahedron will be collapsed. If one edge is collapsed we obtain
+ * a triangle.
+ * @return a Triangle object or NULL
+ */
+ virtual Element* reviseElement() const;
+
protected:
/// Constructor without nodes (for use of derived classes)
Tet(unsigned value = 0);
@@ -87,8 +98,13 @@ protected:
/// Calculates the volume of a tetrahedron via the determinant of the matrix given by its four points.
double computeVolume();
- /// Return a specific edge node.
- inline Node* getEdgeNode(unsigned edge_id, unsigned node_id) const { return _nodes[_edge_nodes[edge_id][node_id]]; };
+ /**
+ * Return a specific edge node, see @ref TetrahedronNodeAndEdgeNumbering for numbering
+ * @param edge_id the id/number of the edge, have to be an integer value in the interval [0,6]
+ * @param node_id the id of the node within the edge (either 0 or 1)
+ * @return a pointer to the internal Node
+ */
+ inline Node const* getEdgeNode(unsigned edge_id, unsigned node_id) const { return _nodes[_edge_nodes[edge_id][node_id]]; };
static const unsigned _face_nodes[4][3];
static const unsigned _edge_nodes[6][2];
diff --git a/MeshLib/Elements/Tri.cpp b/MeshLib/Elements/Tri.cpp
index 3778b7c0e160609ddc2577364d04f7c443b731b0..653343ca2a3c43d63876252fa971f87d7262516a 100644
--- a/MeshLib/Elements/Tri.cpp
+++ b/MeshLib/Elements/Tri.cpp
@@ -11,6 +11,7 @@
*/
#include "Tri.h"
+#include "Edge.h"
#include "Node.h"
#include "MathTools.h"
@@ -76,5 +77,19 @@ double Tri::computeArea()
return MathLib::calcTriangleArea(_nodes[0]->getCoords(), _nodes[1]->getCoords(), _nodes[2]->getCoords());
}
+Element* Tri::reviseElement() const
+{
+ // try to create an edge
+ if (_nodes[0] == _nodes[1] || _nodes[1] == _nodes[2]) {
+ return new Edge(_nodes[0], _nodes[2], _value);
+ }
+
+ if (_nodes[0] == _nodes[2]) {
+ return new Edge(_nodes[0], _nodes[1], _value);
+ }
+
+ return NULL;
+}
+
}
diff --git a/MeshLib/Elements/Tri.h b/MeshLib/Elements/Tri.h
index 3440323edb945e5c9068379babfc1a6fda810f18..d6ebcfce0538de9b3b5b02368fc9274ab536382f 100644
--- a/MeshLib/Elements/Tri.h
+++ b/MeshLib/Elements/Tri.h
@@ -18,18 +18,19 @@
namespace MeshLib {
/**
- * A 2d Triangle Element.
+ * This class represents a 2d triangle element. The following sketch shows the node and edge numbering.
+ * @anchor TriNodeAndEdgeNumbering
* @code
*
- * Tri: 2
+ * 2
* o
* / \
* / \
- * / \
+ * 2/ \1
* / \
* / \
- * o-----------o
- * 0 1
+ * 0-----------1
+ * 0
*
* @endcode
*/
@@ -69,12 +70,22 @@ public:
*/
virtual Element* clone() const;
+ /**
+ * This method should be called after at least two nodes of the triangle
+ * element are collapsed. As a consequence of the node collapsing an edge
+ * of the triangle will be collapsed. If one of the edges is collapsed we
+ * obtain an edge. In this case the method will create the appropriate
+ * object of class Edge.
+ * @return an Edge object or NULL
+ */
+ virtual Element* reviseElement() const;
+
protected:
/// Calculates the area of the triangle by returning half of the area of the corresponding parallelogram.
double computeArea();
/// Return a specific edge node.
- inline Node* getEdgeNode(unsigned edge_id, unsigned node_id) const { return _nodes[_edge_nodes[edge_id][node_id]]; };
+ inline Node const* getEdgeNode(unsigned edge_id, unsigned node_id) const { return _nodes[_edge_nodes[edge_id][node_id]]; };
static const unsigned _edge_nodes[3][2];
diff --git a/MeshLib/MeshCoarsener.cpp b/MeshLib/MeshCoarsener.cpp
index 8abd103d3f7b1320c5a9649247b8a578d02460e0..a739714207ffc0d39a6bcf33d2dcde412ad4ec45 100644
--- a/MeshLib/MeshCoarsener.cpp
+++ b/MeshLib/MeshCoarsener.cpp
@@ -9,6 +9,9 @@
* Created on Aug 3, 2012 by Thomas Fischer
*/
+// BaseLib/logog
+#include "logog.hpp"
+
#include "MeshCoarsener.h"
// GeoLib
@@ -56,7 +59,7 @@ Mesh* MeshCoarsener::operator()(double min_distance)
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());
+ const size_t node_id_k(node->getID());
grid->getVecsOfGridCellsIntersectingCube(node->getCoords(), min_distance, node_vecs_intersecting_cube);
const size_t n_vecs (node_vecs_intersecting_cube.size());
@@ -64,10 +67,15 @@ Mesh* MeshCoarsener::operator()(double min_distance)
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) {
+ Node const*const test_node((*node_vec)[j]);
+ const size_t test_node_id (test_node->getID());
+ if (node_id_k < test_node_id) {
+ if (MathLib::sqrDist(node->getCoords(), test_node->getCoords()) < sqr_min_distance) {
// two nodes are very close to each other
- id_map[j] = k;
+ id_map[test_node_id] = node_id_k;
+#ifndef NDEBUG
+ INFO ("distance of nodes with ids %d and %d is %f", node_id_k, test_node_id, sqrt(MathLib::sqrDist(node->getCoords(), test_node->getCoords())));
+#endif
}
}
}
@@ -86,10 +94,12 @@ Mesh* MeshCoarsener::operator()(double min_distance)
}
// delete unused nodes
- for (size_t k(0); k < n_nodes; k++) {
- if (id_map[k] != k) {
+ for (size_t k(0), cnt(0); k < n_nodes; k++) {
+ if (id_map[k] != cnt) {
delete nodes[k];
nodes[k] = NULL;
+ } else {
+ cnt++;
}
}
@@ -102,26 +112,70 @@ Mesh* MeshCoarsener::operator()(double min_distance)
}
}
+ // reset mesh node ids
+ const size_t new_n_nodes(nodes.size());
+ for (size_t k(0); k < new_n_nodes; k++) {
+ nodes[k]->setID(k);
+ }
+
// 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());
+ std::vector<size_t> mapped_node_ids_of_element;
+ for (size_t k(0), cnt(0); k < n_elements; k++) {
+ Element const*const kth_orig_elem(orig_elements[k]);
+ const size_t n_nodes_element (kth_orig_elem->getNNodes());
+ mapped_node_ids_of_element.clear();
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());
+ const size_t orig_node_id (kth_orig_elem->getNode(i)->getID());
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;
+ } else {
+ mapped_node_ids_of_element.push_back(id_map[it->second]);
+ }
+ }
+
+ // check if nodes of the element are collapsed
+ bool not_collapsed (true);
+ for (size_t i(0); i<n_nodes_element-1 && not_collapsed; i++) {
+ const size_t id_i(mapped_node_ids_of_element[i]);
+ for (size_t j(i+1); j<n_nodes_element && not_collapsed; j++) {
+ if (id_i == mapped_node_ids_of_element[j]) {
+ not_collapsed = false;
+ }
}
- orig_node_pos = it->second;
- elements[k]->setNode(i, nodes[id_map[orig_node_pos]]);
+ }
+ if (! not_collapsed) {
+ Element* elem (kth_orig_elem->clone());
+ if (elem != NULL) {
+ for (size_t i(0); i<n_nodes_element; i++) {
+ elem->setNode(i, nodes[mapped_node_ids_of_element[i]]);
+ }
+ Element* revised_elem(elem->reviseElement());
+ elements[cnt] = revised_elem;
+ delete elem;
+ cnt++;
+ }
+ } else {
+ elements[cnt] = kth_orig_elem->clone();
+ for (size_t i(0); i<n_nodes_element; i++) {
+ elements[cnt]->setNode(i, nodes[mapped_node_ids_of_element[i]]);
+ }
+ cnt++;
+ }
+ }
+
+ for (std::vector<Element*>::iterator it(elements.begin()); it != elements.end(); ) {
+ if (*it == NULL) {
+ it = elements.erase(it);
+ } else {
+ it++;
}
}
- return new Mesh ("test", nodes, elements);
+ return new Mesh (_orig_mesh->getName() + "Collapsed", nodes, elements);
}
} // end namespace MeshLib
diff --git a/MeshLib/Node.h b/MeshLib/Node.h
index 703afbb8d2c3ed2cc8ef9b9fb3e7ce712c3444f2..740c05bb58959020ebbab69226da01ec9b900759 100644
--- a/MeshLib/Node.h
+++ b/MeshLib/Node.h
@@ -29,8 +29,9 @@ class Element;
*/
class Node : public GeoLib::PointWithID
{
- /* friend functions: */
+ /* friend classes: */
friend class Mesh;//void Mesh::setElementInformationForNodes();
+ friend class MeshCoarsener;
//friend void Mesh::addElement(Element*);
diff --git a/SimpleTests/MeshTests/CollapseMeshNodes.cpp b/SimpleTests/MeshTests/CollapseMeshNodes.cpp
index 9e901b5860db5142cbbbfe36b2b1b3a71ed075a0..ee44736e2088a4aa288366480f328231f6cb9542 100644
--- a/SimpleTests/MeshTests/CollapseMeshNodes.cpp
+++ b/SimpleTests/MeshTests/CollapseMeshNodes.cpp
@@ -24,22 +24,42 @@
#include "Legacy/MeshIO.h"
#include "MeshCoarsener.h"
+
+/**
+ * new formatter for logog
+ */
+class FormatterCustom : public logog::FormatterGCC
+{
+ virtual TOPIC_FLAGS GetTopicFlags( const logog::Topic &topic )
+ {
+ return ( Formatter::GetTopicFlags( topic ) &
+ ~( TOPIC_FILE_NAME_FLAG | TOPIC_LINE_NUMBER_FLAG ));
+ }
+};
+
int main(int argc, char *argv[])
{
LOGOG_INITIALIZE();
- logog::Cout* logogCout = new logog::Cout;
+ FormatterCustom *custom_format (new FormatterCustom);
+ logog::Cout *logogCout(new logog::Cout);
+ logogCout->SetFormatter(*custom_format);
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");
-
+ TCLAP::ValueArg<std::string> input_mesh_arg("m","mesh","input mesh file name",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 );
+ TCLAP::ValueArg<std::string> output_mesh_arg("","out-mesh","mesh file name for output",false,"","string");
+ cmd.add( output_mesh_arg );
+
+ TCLAP::ValueArg<double> distance_arg("d","collapse-distance","maximal distance two nodes are collapsed",false,0.01,"for example you can set this parameter to 10^{-6} times maximal area length");
+ cmd.add( distance_arg );
+
cmd.parse( argc, argv );
std::string fname (input_mesh_arg.getValue());
@@ -55,11 +75,9 @@ int main(int argc, char *argv[])
#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());
}
@@ -70,7 +88,8 @@ int main(int argc, char *argv[])
run_time.start();
#endif
MeshLib::MeshCoarsener mesh_coarsener(mesh);
- MeshLib::Mesh *collapsed_mesh(mesh_coarsener (10));
+ MeshLib::Mesh *collapsed_mesh(mesh_coarsener (distance_arg.getValue()));
+
#ifndef WIN32
run_time.stop();
unsigned long mem_with_meshgrid (mem_watch.getVirtMemUsage());
@@ -79,13 +98,19 @@ int main(int argc, char *argv[])
#endif
mesh_io.setMesh(collapsed_mesh);
- std::string out_fname("/home/fischeth/workspace/OGS-6/Build/CollapsedMesh.msh");
+ std::string out_fname (output_mesh_arg.getValue());
+ if (out_fname.empty()) {
+ 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 custom_format;
delete logogCout;
LOGOG_SHUTDOWN();
+
+ return 0;
}