diff --git a/.travis.yml b/.travis.yml
index 6cfe033aed4aa10e83ddd09d53e73875f9d58908..c54b27a0e255b5fbfbce4c0adf52f32ebe9e0685 100644
--- a/.travis.yml
+++ b/.travis.yml
@@ -40,6 +40,8 @@ before_install:
script:
- "pwd & mkdir build && cd build && cmake $CMAKE_ARGS .. && cmake .. && make"
- make test
+ # PetSc
+ - if [[ "$CASE" == "CLI_PETSC" ]]; then make test_mpi; fi
notifications:
email:
diff --git a/FileIO/XmlIO/Boost/BoostXmlGmlInterface.cpp b/FileIO/XmlIO/Boost/BoostXmlGmlInterface.cpp
index 6a591087c90c2813337132a29e0c3dfc8063a5ce..54e47f514b8460a1603e3b566d7b30929e110e14 100644
--- a/FileIO/XmlIO/Boost/BoostXmlGmlInterface.cpp
+++ b/FileIO/XmlIO/Boost/BoostXmlGmlInterface.cpp
@@ -28,8 +28,8 @@
namespace FileIO
{
-BoostXmlGmlInterface::BoostXmlGmlInterface(ProjectData & project_data) :
- _project_data(project_data)
+BoostXmlGmlInterface::BoostXmlGmlInterface(GeoLib::GEOObjects& geo_objs) :
+ _geo_objects(geo_objs)
{}
bool BoostXmlGmlInterface::readFile(const std::string &fname)
@@ -50,7 +50,7 @@ bool BoostXmlGmlInterface::readFile(const std::string &fname)
std::map<std::string, std::size_t>* ply_names = new std::map<std::string, std::size_t>;
std::map<std::string, std::size_t>* sfc_names = new std::map<std::string, std::size_t>;
- GeoLib::GEOObjects* geo_objects (_project_data.getGEOObjects());
+ GeoLib::GEOObjects* geo_objects (&_geo_objects);
// build DOM tree
diff --git a/FileIO/XmlIO/Boost/BoostXmlGmlInterface.h b/FileIO/XmlIO/Boost/BoostXmlGmlInterface.h
index 197a3f54500a3ea7881049bbeba6826fb3f27270..f79a873982ce7af0da0fbcfbff1f581fd6082427 100644
--- a/FileIO/XmlIO/Boost/BoostXmlGmlInterface.h
+++ b/FileIO/XmlIO/Boost/BoostXmlGmlInterface.h
@@ -25,7 +25,6 @@
#include "../XMLInterface.h"
-class ProjectData;
class Point;
class Polyline;
class Surface;
@@ -36,7 +35,7 @@ namespace FileIO
class BoostXmlGmlInterface : public XMLInterface
{
public:
- BoostXmlGmlInterface(ProjectData & project);
+ BoostXmlGmlInterface(GeoLib::GEOObjects& geo_objs);
virtual ~BoostXmlGmlInterface() {}
/// Reads an xml-file containing OGS geometry
@@ -70,7 +69,7 @@ private:
bool isGmlFile( boost::property_tree::ptree const& root) const;
std::map<std::size_t, std::size_t> _idx_map;
- ProjectData & _project_data;
+ GeoLib::GEOObjects &_geo_objects;
};
} // end namespace FileIO
diff --git a/GeoLib/AnalyticalGeometry.cpp b/GeoLib/AnalyticalGeometry.cpp
index 37771bc670bbc9f4d9a0b1dc7d929204b01cdd22..362b8451d8cd9ea3eee7ea96e578c8ad479c1148 100644
--- a/GeoLib/AnalyticalGeometry.cpp
+++ b/GeoLib/AnalyticalGeometry.cpp
@@ -103,6 +103,19 @@ bool lineSegmentIntersect(
if (!isCoplanar(a, b, c, d))
return false;
+ // handle special cases here to avoid computing intersection numerical
+ if (MathLib::sqrDist(a, c) < std::numeric_limits<double>::epsilon() ||
+ MathLib::sqrDist(a, d) < std::numeric_limits<double>::epsilon()) {
+ s = a;
+ return true;
+ }
+ if (MathLib::sqrDist(b, c) < std::numeric_limits<double>::epsilon() ||
+ MathLib::sqrDist(b, d) < std::numeric_limits<double>::epsilon()) {
+ s = b;
+ return true;
+ }
+
+ // general case
MathLib::Vector3 const v(a, b);
MathLib::Vector3 const w(c, d);
MathLib::Vector3 const qp(a, c);
@@ -155,7 +168,7 @@ bool lineSegmentIntersect(
return false;
}
-bool lineSegmentsIntersect(const GeoLib::Polyline* ply,
+bool lineSegmentsIntersect(const GeoLib::Polyline* ply,
size_t &idx0,
size_t &idx1,
GeoLib::Point& intersection_pnt)
@@ -184,70 +197,59 @@ bool lineSegmentsIntersect(const GeoLib::Polyline* ply,
return false;
}
-static
-bool isPointInTriangle(const double p[3], const double a[3], const double b[3], const double c[3])
-{
- // criterion: p-b = u0 * (b - a) + u1 * (b - c); 0 <= u0, u1 <= 1, u0+u1 <= 1
- MathLib::DenseMatrix<double> mat(2, 2);
- mat(0, 0) = a[0] - b[0];
- mat(0, 1) = c[0] - b[0];
- mat(1, 0) = a[1] - b[1];
- mat(1, 1) = c[1] - b[1];
- double rhs[2] = { p[0] - b[0], p[1] - b[1] };
-
- MathLib::GaussAlgorithm<MathLib::DenseMatrix<double>, double*> gauss(mat);
- gauss.solve(rhs);
-
- if (0 <= rhs[0] && rhs[0] <= 1 && 0 <= rhs[1] && rhs[1] <= 1 && rhs[0] + rhs[1] <= 1)
- return true;
- return false;
-}
-
bool isPointInTriangle(const GeoLib::Point* p, const GeoLib::Point* a, const GeoLib::Point* b,
const GeoLib::Point* c)
{
- return isPointInTriangle(p->getCoords(), a->getCoords(), b->getCoords(), c->getCoords());
+ return isPointInTriangle(*p, *a, *b, *c);
}
-static
-double getOrientedTriArea(GeoLib::Point const& a, GeoLib::Point const& b, GeoLib::Point const& c)
+bool isPointInTriangle(GeoLib::Point const& q,
+ GeoLib::Point const& a,
+ GeoLib::Point const& b,
+ GeoLib::Point const& c,
+ double eps)
{
- const MathLib::Vector3 u(a,c);
- const MathLib::Vector3 v(a,b);
- const MathLib::Vector3 w(MathLib::crossProduct(u, v));
- return 0.5 * sqrt(MathLib::scalarProduct(w, w));
-}
+ MathLib::Vector3 const v(a, b);
+ MathLib::Vector3 const w(a, c);
-bool isPointInTriangle(GeoLib::Point const& p, GeoLib::Point const& a, GeoLib::Point const& b,
- GeoLib::Point const& c, double eps)
-{
- const unsigned dim(3);
- MathLib::DenseMatrix<double> m(dim, dim);
- for (unsigned i(0); i < dim; i++)
- m(i, 0) = b[i] - a[i];
- for (unsigned i(0); i < dim; i++)
- m(i, 1) = c[i] - a[i];
- for (unsigned i(0); i < dim; i++)
- m(i, 2) = p[i] - a[i];
-
- // point p is in the same plane as the triangle if and only if
- // the following determinate of the 3x3 matrix equals zero (up to an eps)
- double det3x3(m(0, 0) * (m(1, 1) * m(2, 2) - m(2, 1) * m(1, 2))
- - m(1, 0) * (m(2, 1) * m(0, 2) - m(0, 1) * m(2, 2))
- + m(2, 0) * (m(0, 1) * m(1, 2) - m(1, 1) * m(0, 2)));
- if (fabs(det3x3) > eps)
- return false;
+ MathLib::DenseMatrix<double> mat (2,2);
+ mat(0,0) = v.getSqrLength();
+ mat(0,1) = v[0]*w[0] + v[1]*w[1] + v[2]*w[2];
+ mat(1,0) = mat(0,1);
+ mat(1,1) = w.getSqrLength();
+ double y[2] = {
+ v[0] * (q[0] - a[0]) + v[1] * (q[1] - a[1]) + v[2] * (q[2] - a[2]),
+ w[0] * (q[0] - a[0]) + w[1] * (q[1] - a[1]) + w[2] * (q[2] - a[2])
+ };
- double total_area(getOrientedTriArea(a, b, c));
- double abp_area(getOrientedTriArea(a, b, p));
- double bcp_area(getOrientedTriArea(b, c, p));
- double cap_area(getOrientedTriArea(c, a, p));
+ MathLib::GaussAlgorithm<MathLib::DenseMatrix<double>, double*> gauss(mat);
+ gauss.solve(y);
+
+ const double lower (std::numeric_limits<float>::epsilon());
+ const double upper (1 + lower);
+
+ if (-lower <= y[0] && y[0] <= upper && -lower <= y[1] && y[1] <= upper && y[0] + y[1] <=
+ upper) {
+ GeoLib::Point const q_projected(
+ a[0] + y[0] * v[0] + y[1] * w[0],
+ a[1] + y[0] * v[1] + y[1] * w[1],
+ a[2] + y[0] * v[2] + y[1] * w[2]
+ );
+ if (MathLib::sqrDist(q, q_projected) < eps)
+ return true;
+ }
- if (fabs(abp_area + bcp_area + cap_area - total_area) < eps)
- return true;
return false;
}
+double calcTriangleArea(GeoLib::Point const& a, GeoLib::Point const& b, GeoLib::Point const& c)
+{
+ MathLib::Vector3 const u(a,c);
+ MathLib::Vector3 const v(a,b);
+ MathLib::Vector3 const w(MathLib::crossProduct(u, v));
+ return 0.5 * w.getLength();
+}
+
// NewellPlane from book Real-Time Collision detection p. 494
void getNewellPlane(const std::vector<GeoLib::Point*>& pnts, MathLib::Vector3 &plane_normal, double& d)
{
@@ -435,4 +437,26 @@ bool isCoplanar(const GeoLib::Point& a, const GeoLib::Point& b, const GeoLib::Po
return (sqr_scalar_triple/normalisation_factor < 1e-11);
}
+void computeAndInsertAllIntersectionPoints(GeoLib::PointVec &pnt_vec,
+ std::vector<GeoLib::Polyline*> & plys)
+{
+ for (auto it0(plys.begin()); it0 != plys.end(); ++it0) {
+ auto it1(it0);
+ ++it1;
+ for (; it1 != plys.end(); ++it1) {
+ GeoLib::Point s;
+ for (std::size_t i(0); i<(*it0)->getNumberOfPoints()-1; i++) {
+ for (std::size_t j(0); j<(*it1)->getNumberOfPoints()-1; j++) {
+ if (lineSegmentIntersect(*(*it0)->getPoint(i), *(*it0)->getPoint(i+1),
+ *(*it1)->getPoint(j), *(*it1)->getPoint(j+1), s)) {
+ std::size_t const id(pnt_vec.push_back(new GeoLib::Point(s)));
+ (*it0)->insertPoint(i+1, id);
+ (*it1)->insertPoint(j+1, id);
+ }
+ }
+ }
+ }
+ }
+}
+
} // end namespace GeoLib
diff --git a/GeoLib/AnalyticalGeometry.h b/GeoLib/AnalyticalGeometry.h
index 9aed682d823bf508076880a638387cc83735f726..9a705b5d1e2828dd27a2c94107a0a7b99696bc3a 100644
--- a/GeoLib/AnalyticalGeometry.h
+++ b/GeoLib/AnalyticalGeometry.h
@@ -17,6 +17,7 @@
// GeoLib
#include "Triangle.h"
+#include "PointVec.h"
// MathLib
#include "LinAlg/Dense/DenseMatrix.h"
@@ -98,12 +99,29 @@ void rotatePointsToXY(std::vector<GeoLib::Point*> &pnts);
*/
void rotatePointsToXZ(std::vector<GeoLib::Point*> &pnts);
+/**
+ * Calculates the area of the triangle defined by its edge nodes a, b and c..
+ * The formula is \f$A= \frac{1}{2} \cdot |u \times v|\f$, i.e. half of the area of the
+ * parallelogram specified by the vectors\f$u=b-a\f$ and \f$v=c-a\f.
+ */
+double calcTriangleArea(GeoLib::Point const& a, GeoLib::Point const& b, GeoLib::Point const& c);
+
bool isPointInTriangle (const GeoLib::Point* p,
const GeoLib::Point* a, const GeoLib::Point* b, const GeoLib::Point* c);
+/**
+ * Tests if the given point p is within the triangle, defined by its edge nodes a, b and c.
+ * Using the eps it is possible to test a 'epsilon' neighbourhood around the triangle.
+ * @param p test point
+ * @param a edge node of triangle
+ * @param b edge node of triangle
+ * @param c edge node of triangle
+ * @param eps size of neighbourhood (orthogonal distance to the plane spaned by triangle)
+ * @return true if the test point p is within the 'epsilon'-neighbourhood of the triangle
+ */
bool isPointInTriangle(GeoLib::Point const& p,
GeoLib::Point const& a, GeoLib::Point const& b, GeoLib::Point const& c,
- double eps = std::numeric_limits<double>::epsilon());
+ double eps = std::numeric_limits<float>::epsilon());
/**
* test for intersections of the line segments of the Polyline
@@ -166,6 +184,16 @@ double scalarTriple(MathLib::Vector3 const& u, MathLib::Vector3 const& v, MathLi
const GeoLib::Point& c, const GeoLib::Point& d);
+/**
+ * Method first computes the intersection points of line segements of GeoLib::Polyline objects
+ * (@see computeIntersectionPoints()) and pushes each intersection point in the GeoLib::PointVec
+ * pnt_vec. For each intersection an id is returned. This id is used to split the two
+ * intersecting straight line segments in four straight line segments.
+ */
+void computeAndInsertAllIntersectionPoints(
+ GeoLib::PointVec &pnt_vec,
+ std::vector<GeoLib::Polyline*> & plys);
+
} // end namespace GeoLib
#endif /* ANALYTICAL_GEOMETRY_H_ */
diff --git a/GeoLib/Polygon.cpp b/GeoLib/Polygon.cpp
index 93c4c394ea84f9df49d1fe0558c4e69a9bd3c3e5..fc3ebfe3e76581ea46b5738c84cfd0ec6e8eebe6 100644
--- a/GeoLib/Polygon.cpp
+++ b/GeoLib/Polygon.cpp
@@ -113,18 +113,84 @@ bool Polygon::isPntInPolygon(double x, double y, double z) const
return isPntInPolygon (pnt);
}
-bool Polygon::isPolylineInPolygon(const Polyline& ply) const
+std::vector<GeoLib::Point> Polygon::getAllIntersectionPoints(
+ GeoLib::Point const& a, GeoLib::Point const& b) const
{
- std::size_t ply_size (ply.getNumberOfPoints()), cnt (0);
- for (std::size_t k(0); k < ply_size; k++) {
- if (isPntInPolygon (*(ply.getPoint(k)))) {
- cnt++;
+ std::vector<GeoLib::Point> intersections;
+ const std::size_t n_segments(getNumberOfPoints() - 1);
+ GeoLib::Point s;
+ for (std::size_t k(0); k < n_segments; k++) {
+ if (GeoLib::lineSegmentIntersect(*(getPoint(k)), *(getPoint(k+1)), a, b, s)) {
+ intersections.push_back(s);
}
}
- if (cnt == ply_size)
- return true;
- return false;
+ return intersections;
+}
+
+bool Polygon::containsSegment(GeoLib::Point const& a, GeoLib::Point const& b) const
+{
+ std::vector<GeoLib::Point> s(getAllIntersectionPoints(a, b));
+
+ // no intersections -> check if at least one point of segment is in polygon
+ if (s.empty()) {
+ return (isPntInPolygon(a));
+ }
+
+ const double tol(std::numeric_limits<float>::epsilon());
+
+ // one intersection, intersection in line segment end point
+ if (s.size() == 1) {
+ const double sqr_dist_as(MathLib::sqrDist(a,s[0]));
+ if (sqr_dist_as < tol) {
+ return (isPntInPolygon(b));
+ }
+
+ const double sqr_dist_bs(MathLib::sqrDist(b,s[0]));
+ if (sqr_dist_bs < tol) {
+ return (isPntInPolygon(a));
+ }
+ }
+
+ // Sorting the intersection with respect to the distance to the point a.
+ // This induces a partition of the line segment into sub segments.
+ std::sort(s.begin(), s.end(),
+ [&a] (GeoLib::Point const& p0, GeoLib::Point const& p1) {
+ return MathLib::sqrDist(a, p0) <= MathLib::sqrDist(a, p1);
+ }
+ );
+
+ // remove sub segments with almost zero length
+ for (std::size_t k(0); k<s.size()-1; ) {
+ if (MathLib::sqrDist(s[k], s[k+1]) < tol) {
+ s.erase(s.begin()+k+1);
+ } else {
+ k++;
+ }
+ }
+
+ // Check if all sub segments are within the polygon.
+ if (!isPntInPolygon(GeoLib::Point(0.5*(a[0]+s[0][0]), 0.5*(a[1]+s[0][1]), 0.5*(a[2]+s[0][2]))))
+ return false;
+ const std::size_t n_sub_segs(s.size()-1);
+ for (std::size_t k(0); k<n_sub_segs; k++) {
+ if (!isPntInPolygon(GeoLib::Point(0.5*(s[k][0]+s[k+1][0]), 0.5*(s[k][1]+s[k+1][1]), 0.5*(s[k][2]+s[k+1][2]))))
+ return false;
+ }
+ if (!isPntInPolygon(GeoLib::Point(0.5*(s[0][0]+b[0]), 0.5*(s[0][1]+b[1]), 0.5*(s[0][2]+b[2]))))
+ return false;
+ return true;
+}
+
+bool Polygon::isPolylineInPolygon(const Polyline& ply) const
+{
+ std::size_t const n_segments(ply.getNumberOfPoints()-1);
+ for (std::size_t k(0); k < n_segments; k++) {
+ if (!containsSegment(*ply.getPoint(k), *ply.getPoint(k+1))) {
+ return false;
+ }
+ }
+ return true;
}
bool Polygon::isPartOfPolylineInPolygon(const Polyline& ply) const
diff --git a/GeoLib/Polygon.h b/GeoLib/Polygon.h
index 667d935d8602cc178cd0f126a58da57814e89224..d66dd5998bad4111be2e6fef6bfb8e7016416512 100644
--- a/GeoLib/Polygon.h
+++ b/GeoLib/Polygon.h
@@ -39,7 +39,7 @@ enum class EdgeType
};
/**
- *
+ * A polygon is a (closed) polyline. Thus class Polygon is derived from class Polyline.
*/
class Polygon : public Polyline
{
@@ -71,6 +71,16 @@ public:
* @return if point is inside the polygon true, else false
*/
bool isPntInPolygon (double x, double y, double z) const;
+
+ /**
+ * Checks if the straight line segment given by its end points a and b
+ * is contained within the polygon.
+ * @param a the first end point of the straight line segment
+ * @param b the second end point of the straight line segment
+ * @return true if the straight line segment is within the polygon, else false
+ */
+ bool containsSegment(GeoLib::Point const& a, GeoLib::Point const& b) const;
+
/**
* Method checks if all points of the polyline ply are inside of the polygon.
* @param ply the polyline that should be checked
@@ -104,6 +114,16 @@ public:
friend bool operator==(Polygon const& lhs, Polygon const& rhs);
private:
+ /**
+ * computes all intersection points of the straight line segment and the polyline boundary
+ * @param a end point of line segment
+ * @param b end point of line segment
+ * @return a vector of tuples, where a tuple contains the intersection point and
+ * the intersected segment number
+ */
+ std::vector<GeoLib::Point> getAllIntersectionPoints(
+ GeoLib::Point const& a, GeoLib::Point const& b) const;
+
/**
* from book: Computational Geometry and Computer Graphics in C++, page 119
* get the type of edge with respect to the given point (2d method!)
diff --git a/GeoLib/PolylineWithSegmentMarker.h b/GeoLib/PolylineWithSegmentMarker.h
index 0e7cc7a905d800bf5bc9affe882cce5e7eee243d..f904bf032af4b28d19b1f8f7c714427d6131b2d0 100644
--- a/GeoLib/PolylineWithSegmentMarker.h
+++ b/GeoLib/PolylineWithSegmentMarker.h
@@ -19,6 +19,10 @@
namespace GeoLib {
+/**
+ * This is a polyline with the possibility to mark some segments. Thus class
+ * PolylineWithSegmentMarker is derived from class Polyline.
+ */
class PolylineWithSegmentMarker: public GeoLib::Polyline {
public:
PolylineWithSegmentMarker(GeoLib::Polyline const& polyline);
diff --git a/GeoLib/SimplePolygonTree.cpp b/GeoLib/SimplePolygonTree.cpp
index d0e5835226b254e4d890cc76f7dc0de48ac55e7f..6c09108b1e2aeb47b7461cbddc37e1734554b768 100644
--- a/GeoLib/SimplePolygonTree.cpp
+++ b/GeoLib/SimplePolygonTree.cpp
@@ -22,7 +22,6 @@ SimplePolygonTree::SimplePolygonTree(Polygon * polygon, SimplePolygonTree * pare
SimplePolygonTree::~SimplePolygonTree()
{
- delete _node_polygon;
for (std::list<SimplePolygonTree*>::const_iterator it (_childs.begin());
it != _childs.end(); ++it) {
delete *it;
@@ -31,20 +30,7 @@ SimplePolygonTree::~SimplePolygonTree()
bool SimplePolygonTree::isPolygonInside (const SimplePolygonTree* polygon_hierarchy) const
{
- const Polygon* polygon (polygon_hierarchy->getPolygon());
- // check *all* points of polygon
- size_t n_pnts_polygon (polygon->getNumberOfPoints() - 1), cnt(0);
- for (size_t k(0); k < n_pnts_polygon && cnt == k; k++) {
- if (_node_polygon->isPntInPolygon (*(polygon->getPoint(k)))) {
- cnt++;
- }
- }
-
- // all points of the given polygon are contained in the
- if (cnt == n_pnts_polygon)
- return true;
- else
- return false;
+ return _node_polygon->isPolylineInPolygon(*(polygon_hierarchy->getPolygon()));
}
void SimplePolygonTree::insertSimplePolygonTree (SimplePolygonTree* polygon_hierarchy)
@@ -53,15 +39,7 @@ void SimplePolygonTree::insertSimplePolygonTree (SimplePolygonTree* polygon_hier
bool nfound (true);
for (std::list<SimplePolygonTree*>::const_iterator it (_childs.begin());
it != _childs.end() && nfound; ++it) {
- // check all points of polygon
- size_t n_pnts_polygon (polygon->getNumberOfPoints()), cnt(0);
- for (size_t k(0); k < n_pnts_polygon && cnt == k; k++) {
- if (((*it)->getPolygon())->isPntInPolygon (*(polygon->getPoint(k)))) {
- cnt++;
- }
- }
- // all points of the given polygon are contained in the current polygon
- if (cnt == n_pnts_polygon) {
+ if (((*it)->getPolygon())->isPolylineInPolygon (*(polygon))) {
(*it)->insertSimplePolygonTree (polygon_hierarchy);
nfound = false;
}
diff --git a/GeoLib/SimplePolygonTree.h b/GeoLib/SimplePolygonTree.h
index b5212bb50b150167cf5475c9b645e1350c5a240f..adbaf2a157a5fc913e0824823454dd3d0614fba7 100644
--- a/GeoLib/SimplePolygonTree.h
+++ b/GeoLib/SimplePolygonTree.h
@@ -24,16 +24,28 @@ namespace GeoLib
{
/**
* \brief This class computes and stores the topological relations between
- * polygons. Every node of the SimplePolygonTree represents a polygon.
- *
+ * polygons. Every node of the SimplePolygonTree represents a polygon. A
+ * child node c of a parent node p mean that the polygon represented by c
+ * is contained in the polygon represented by p.
*/
class SimplePolygonTree
{
public:
+ /** Creates a node of a tree containing a simple polygon.
+ * @param polygon the polygon represented by this tree node
+ * @param parent pointer to the parent node within the tree or nullptr
+ * (if SimplePolygonTree node is the root node of the tree)
+ */
SimplePolygonTree(Polygon* polygon, SimplePolygonTree* parent);
+ /** Destructor: Attention: does not destroy the polygon! */
virtual ~SimplePolygonTree();
-
+ /** Checks if the polygon represented by the given polygon tree node
+ * is inside this node polygon.
+ */
bool isPolygonInside (const SimplePolygonTree* polygon_tree) const;
+ /** Either insert the given SimplePolygonTree in one of the existing
+ * childs or as a new child.
+ */
void insertSimplePolygonTree (SimplePolygonTree* polygon_tree);
/**
@@ -42,6 +54,7 @@ public:
*/
const Polygon* getPolygon () const;
+ /** returns the number of childs */
std::size_t getNChilds() const { return _childs.size(); }
protected:
@@ -59,6 +72,7 @@ protected:
* in the _node_polygon
*/
std::list<SimplePolygonTree*> _childs;
+
private:
void setParent(SimplePolygonTree* parent)
{
@@ -72,24 +86,25 @@ private:
template <typename POLYGONTREETYPE>
void createPolygonTrees (std::list<POLYGONTREETYPE*>& list_of_simple_polygon_hierarchies)
{
- typename std::list<POLYGONTREETYPE*>::iterator it0 (list_of_simple_polygon_hierarchies.begin()), it1;
- while (it0 != list_of_simple_polygon_hierarchies.end()) {
- it1 = it0;
- it1++;
+ typedef typename std::list<POLYGONTREETYPE*>::const_iterator CIT;
+ typedef typename std::list<POLYGONTREETYPE*>::iterator IT;
+ for (CIT it0(list_of_simple_polygon_hierarchies.begin());
+ it0 != list_of_simple_polygon_hierarchies.end(); ++it0) {
+ IT it1 = list_of_simple_polygon_hierarchies.begin();
while (it1 != list_of_simple_polygon_hierarchies.end()) {
+ if (it0 == it1) { // don't check same polygons
+ ++it1;
+ // skip test if it1 points to the end after increment
+ if (it1 == list_of_simple_polygon_hierarchies.end())
+ break;
+ }
if ((*it0)->isPolygonInside(*it1)) {
(*it0)->insertSimplePolygonTree(*it1);
it1 = list_of_simple_polygon_hierarchies.erase(it1);
} else {
- if ((*it1)->isPolygonInside(*it0)) {
- (*it1)->insertSimplePolygonTree(*it0);
- it0 = list_of_simple_polygon_hierarchies.erase(it0);
- }
-
- it1++;
+ ++it1;
}
}
- it0++;
}
}
diff --git a/GeoLib/TemplateVec.h b/GeoLib/TemplateVec.h
index c7f1afffb62fdde2f979bb89a7e216af822d7edb..8290f87d3ea168d88a67197f06e6bfc9bcd6e271 100644
--- a/GeoLib/TemplateVec.h
+++ b/GeoLib/TemplateVec.h
@@ -18,6 +18,9 @@
#include <algorithm>
#include <stdexcept>
+#include <map>
+#include <vector>
+#include <string>
namespace GeoLib
{
diff --git a/GeoLib/Triangle.cpp b/GeoLib/Triangle.cpp
index 14ea90d525007cc60466d73dc450c49eba242c91..f915d3da8b4f05e5d2f125f0441c287b87cc5682 100644
--- a/GeoLib/Triangle.cpp
+++ b/GeoLib/Triangle.cpp
@@ -14,6 +14,9 @@
#include "Triangle.h"
+// GeoLib
+#include "AnalyticalGeometry.h"
+
// MathLib
#include "LinAlg/Solvers/GaussAlgorithm.h"
#include "MathTools.h"
@@ -59,106 +62,12 @@ void Triangle::setTriangle (size_t pnt_a, size_t pnt_b, size_t pnt_c)
_longest_edge = sqrt (_longest_edge);
}
-bool Triangle::containsPoint (Point const& pnt) const
+bool Triangle::containsPoint(Point const& q, double eps) const
{
- GeoLib::Point const& a_tmp (*(_pnts[_pnt_ids[0]]));
- GeoLib::Point const& b_tmp (*(_pnts[_pnt_ids[1]]));
- GeoLib::Point const& c_tmp (*(_pnts[_pnt_ids[2]]));
-
- GeoLib::Point s(a_tmp);
- for (size_t k(0); k<3; k++) {
- s[k] += b_tmp[k] + c_tmp[k];
- s[k] /= 3.0;
- }
-
- double eps (1e-2);
- GeoLib::Point const a (a_tmp[0] + eps *(a_tmp[0]-s[0]),
- a_tmp[1] + eps *(a_tmp[1]-s[1]),
- a_tmp[2] + eps *(a_tmp[2]-s[2]));
- GeoLib::Point const b (b_tmp[0] + eps *(b_tmp[0]-s[0]),
- b_tmp[1] + eps *(b_tmp[1]-s[1]),
- b_tmp[2] + eps *(b_tmp[2]-s[2]));
- GeoLib::Point const c (c_tmp[0] + eps *(c_tmp[0]-s[0]),
- c_tmp[1] + eps *(c_tmp[1]-s[1]),
- c_tmp[2] + eps *(c_tmp[2]-s[2]));
-
- const double delta (std::numeric_limits<double>::epsilon());
- const double upper (1+delta);
-
- // check special case where points of triangle have the same x-coordinate
- if (fabs(b[0]-a[0]) <= std::numeric_limits<double>::epsilon() &&
- fabs(c[0]-a[0]) <= std::numeric_limits<double>::epsilon()) {
- // all points of triangle have same x-coordinate
- if (fabs(pnt[0]-a[0]) / _longest_edge <= 1e-3) {
- // criterion: p-a = u0 * (b-a) + u1 * (c-a); 0 <= u0, u1 <= 1, u0+u1 <= 1
- MathLib::DenseMatrix<double> mat (2,2);
- mat(0,0) = b[1] - a[1];
- mat(0,1) = c[1] - a[1];
- mat(1,0) = b[2] - a[2];
- mat(1,1) = c[2] - a[2];
- double y[2] = {pnt[1]-a[1], pnt[2]-a[2]};
-
- MathLib::GaussAlgorithm<MathLib::DenseMatrix<double>, double*> gauss (mat);
- gauss.solve (y);
-
- if (-delta <= y[0] && y[0] <= upper && -delta <= y[1] && y[1] <= upper
- && y[0] + y[1] <= upper) {
- return true;
- } else {
- return false;
- }
- } else {
- return false;
- }
- }
-
- // check special case where points of triangle have the same y-coordinate
- if (fabs(b[1]-a[1]) <= std::numeric_limits<double>::epsilon() &&
- fabs(c[1]-a[1]) <= std::numeric_limits<double>::epsilon()) {
- // all points of triangle have same y-coordinate
- if (fabs(pnt[1]-a[1]) / _longest_edge <= 1e-3) {
- // criterion: p-a = u0 * (b-a) + u1 * (c-a); 0 <= u0, u1 <= 1, u0+u1 <= 1
- MathLib::DenseMatrix<double> mat (2,2);
- mat(0,0) = b[0] - a[0];
- mat(0,1) = c[0] - a[0];
- mat(1,0) = b[2] - a[2];
- mat(1,1) = c[2] - a[2];
- double y[2] = {pnt[0]-a[0], pnt[2]-a[2]};
-
- MathLib::GaussAlgorithm<MathLib::DenseMatrix<double>, double*> gauss (mat);
- gauss.solve (y);
-
- if (-delta <= y[0] && y[0] <= upper && -delta <= y[1] && y[1] <= upper && y[0] + y[1] <= upper) {
- return true;
- } else {
- return false;
- }
- } else {
- return false;
- }
- }
-
- // criterion: p-a = u0 * (b-a) + u1 * (c-a); 0 <= u0, u1 <= 1, u0+u1 <= 1
- MathLib::DenseMatrix<double> mat (2,2);
- mat(0,0) = b[0] - a[0];
- mat(0,1) = c[0] - a[0];
- mat(1,0) = b[1] - a[1];
- mat(1,1) = c[1] - a[1];
- double y[2] = {pnt[0]-a[0], pnt[1]-a[1]};
-
- MathLib::GaussAlgorithm<MathLib::DenseMatrix<double>, double*> gauss (mat);
- gauss.solve (y);
-
- // check if the solution fulfills the third equation
- if (fabs((b[2]-a[2]) * y[0] + (c[2]-a[2]) * y[1] - (pnt[2] - a[2])) < 1e-3) {
- if (-delta <= y[0] && y[0] <= upper && -delta <= y[1] && y[1] <= upper &&
- y[0] + y[1] <= upper) {
- return true;
- }
- return false;
- } else {
- return false;
- }
+ GeoLib::Point const& a(*(_pnts[_pnt_ids[0]]));
+ GeoLib::Point const& b(*(_pnts[_pnt_ids[1]]));
+ GeoLib::Point const& c(*(_pnts[_pnt_ids[2]]));
+ return GeoLib::isPointInTriangle(q, a, b, c, eps);
}
bool Triangle::containsPoint2D (Point const& pnt) const
diff --git a/GeoLib/Triangle.h b/GeoLib/Triangle.h
index 3b0dcb07dfd741cd3256315b885b618519af8d89..caa14f01d0a956856e336e2c17d5ea55fa0adb8c 100644
--- a/GeoLib/Triangle.h
+++ b/GeoLib/Triangle.h
@@ -71,11 +71,12 @@ public:
}
/**
- * checks if point is in triangle
- * @param pnt The input point
+ * Checks if point q is within the triangle, uses GeoLib::isPointInTriangle().
+ * @param q The input point.
+ * @param eps Checks the 'epsilon'-neighbourhood
* @return true, if point is in triangle, else false
*/
- bool containsPoint (Point const& pnt) const;
+ bool containsPoint(Point const& q, double eps = std::numeric_limits<float>::epsilon()) const;
/**
* projects the triangle points to the x-y-plane and
diff --git a/Gui/mainwindow.cpp b/Gui/mainwindow.cpp
index bd829ae459b66fc333e6df4f6d4cdf5bc76def1e..15aba1929e039836854a488b85997f4c5261e870 100644
--- a/Gui/mainwindow.cpp
+++ b/Gui/mainwindow.cpp
@@ -1038,6 +1038,9 @@ void MainWindow::callFileConverter() const
{
QSettings settings;
QString converter_path = settings.value("DataExplorerConverterPath").toString();
+#if _WIN32
+ converter_path = QString("\"").append(converter_path).append("\"");
+#endif
if (!converter_path.isEmpty())
system(converter_path.toAscii());
else
diff --git a/MathLib/LinAlg/Dense/DenseMatrix-impl.h b/MathLib/LinAlg/Dense/DenseMatrix-impl.h
index 1bd2e873c1b61c3587a4d30674343c5669cfbea4..0d871635b057c9adee7f85dc4c63b5127e02d12b 100644
--- a/MathLib/LinAlg/Dense/DenseMatrix-impl.h
+++ b/MathLib/LinAlg/Dense/DenseMatrix-impl.h
@@ -276,6 +276,7 @@ template <typename FP_TYPE, typename IDX_TYPE>
void
DenseMatrix<FP_TYPE, IDX_TYPE>::write (std::ostream &out) const
{
+ out << _n_rows << " " << _n_cols << "\n";
for (IDX_TYPE i = 0; i < _n_rows; i++) {
for (IDX_TYPE j = 0; j < _n_cols; j++) {
out << _data[address(i, j)] << "\t";
diff --git a/MathLib/LinAlg/PETSc/PETScLinearSolver.cpp b/MathLib/LinAlg/PETSc/PETScLinearSolver.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..9118c177db6a24813360927608b51ab8cca5f60f
--- /dev/null
+++ b/MathLib/LinAlg/PETSc/PETScLinearSolver.cpp
@@ -0,0 +1,104 @@
+/*!
+ \file PETScLinearSolver.cpp
+ \brief Definition of class PETScLinearSolver, which defines a solver object
+ based on PETSc routines.
+
+ \author Wenqing Wang
+ \version
+ \date Nov 2011 - Mar 2014
+
+ \copyright
+ Copyright (c) 2013, OpenGeoSys Community (http://www.opengeosys.org)
+ Distributed under a Modified BSD License.
+ See accompanying file LICENSE.txt or
+ http://www.opengeosys.org/project/license
+*/
+
+#include "PETScLinearSolver.h"
+
+namespace MathLib
+{
+PETScLinearSolver::PETScLinearSolver(PETScMatrix &A, const std::string &prefix) : _A(A)
+{
+ KSPCreate(PETSC_COMM_WORLD, &_solver);
+
+ KSPGetPC(_solver, &_pc);
+
+ //
+ if ( !prefix.empty() )
+ {
+ KSPSetOptionsPrefix(_solver, prefix.c_str());
+ }
+
+ KSPSetFromOptions(_solver); // set running time option
+}
+
+bool PETScLinearSolver::solve(const PETScVector &b, PETScVector &x)
+{
+// define TEST_MEM_PETSC
+#ifdef TEST_MEM_PETSC
+ PetscLogDouble mem1, mem2;
+ PetscMemoryGetCurrentUsage(&mem1);
+#endif
+
+ KSPSetOperators(_solver, _A.getRawMatrix(), _A.getRawMatrix(), DIFFERENT_NONZERO_PATTERN);
+
+ KSPSolve(_solver, b.getData(), x.getData());
+
+ KSPConvergedReason reason;
+ KSPGetConvergedReason(_solver, &reason);
+
+ bool converged = true;
+ if(reason > 0)
+ {
+ const char *ksp_type;
+ const char *pc_type;
+ KSPGetType(_solver, &ksp_type);
+ PCGetType(_pc, &pc_type);
+
+ PetscPrintf(PETSC_COMM_WORLD,"\n================================================");
+ PetscPrintf(PETSC_COMM_WORLD, "\nLinear solver %s with %s preconditioner",
+ ksp_type, pc_type);
+
+ PetscInt its;
+ KSPGetIterationNumber(_solver, &its);
+ PetscPrintf(PETSC_COMM_WORLD,"\nConvergence in %d iterations.\n", its);
+ PetscPrintf(PETSC_COMM_WORLD,"\n================================================\n");
+ }
+ else if(reason == KSP_DIVERGED_ITS)
+ {
+ PetscPrintf(PETSC_COMM_WORLD, "\nWarning: maximum number of iterations reached.\n");
+ }
+ else
+ {
+ converged = false;
+ if(reason == KSP_DIVERGED_INDEFINITE_PC)
+ {
+ PetscPrintf(PETSC_COMM_WORLD, "\nDivergence because of indefinite preconditioner,");
+ PetscPrintf(PETSC_COMM_WORLD, "\nTry to run again with -pc_factor_shift_positive_definite option.\n");
+ }
+ else if(reason == KSP_DIVERGED_BREAKDOWN_BICG)
+ {
+ PetscPrintf(PETSC_COMM_WORLD, "\nKSPBICG method was detected so the method could not continue to enlarge the Krylov space.");
+ PetscPrintf(PETSC_COMM_WORLD, "\nTry to run again with another solver.\n");
+ }
+ else if(reason == KSP_DIVERGED_NONSYMMETRIC)
+ {
+ PetscPrintf(PETSC_COMM_WORLD, "\nMatrix or preconditioner is unsymmetric but KSP requires symmetric.\n");
+ }
+ else
+ {
+ PetscPrintf(PETSC_COMM_WORLD, "\nDivergence detected, use command option -ksp_monitor or -log_summary to check the details.\n");
+ }
+ }
+
+#ifdef TEST_MEM_PETSC
+ PetscMemoryGetCurrentUsage(&mem2);
+ PetscPrintf(PETSC_COMM_WORLD, "###Memory usage by solver. Before :%f After:%f Increase:%d\n", mem1, mem2, (int)(mem2 - mem1));
+#endif
+
+ return converged;
+}
+
+} //end of namespace
+
diff --git a/MathLib/LinAlg/PETSc/PETScLinearSolver.h b/MathLib/LinAlg/PETSc/PETScLinearSolver.h
new file mode 100644
index 0000000000000000000000000000000000000000..a4454706ca84163136a914d3c8717647b970ec7f
--- /dev/null
+++ b/MathLib/LinAlg/PETSc/PETScLinearSolver.h
@@ -0,0 +1,80 @@
+/*!
+ \file PETScLinearSolver.h
+ \brief Declaration of class PETScLinearSolver, which defines a solver object
+ based on PETSc routines.
+
+ \author Wenqing Wang
+ \version
+ \date Nov 2011 - Sep 2013
+
+ \copyright
+ Copyright (c) 2013, OpenGeoSys Community (http://www.opengeosys.org)
+ Distributed under a Modified BSD License.
+ See accompanying file LICENSE.txt or
+ http://www.opengeosys.org/project/license
+*/
+
+#ifndef PETSCLINEARSOLVER_H_
+#define PETSCLINEARSOLVER_H_
+
+#include<string>
+
+#include "petscksp.h"
+
+#include "PETScMatrix.h"
+#include "PETScVector.h"
+
+namespace MathLib
+{
+
+/*!
+ A class of linear solver based on PETSc rountines.
+
+ All options for KSP and PC must given in the command line.
+*/
+class PETScLinearSolver
+{
+ public:
+
+ /*!
+ Constructor
+ \param A Matrix, cannot be constant.
+ \param prefix Name used to distinguish the options in the command
+ line for this solver. It can be the name of the PDE
+ that owns an instance of this class.
+ */
+ PETScLinearSolver(PETScMatrix &A, const std::string &prefix="");
+
+ ~PETScLinearSolver()
+ {
+ KSPDestroy(&_solver);
+ }
+
+ /*!
+ Solve a system of equations.
+ \param b The right hand side of the equations.
+ \param x The solutions to be solved.
+ \return true: converged, false: diverged.
+ */
+ bool solve(const PETScVector &b, PETScVector &x);
+
+ /// Get number of iterations.
+ PetscInt getNumberOfIterations() const
+ {
+ PetscInt its = 0;
+ KSPGetIterationNumber(_solver, &its);
+ return its;
+ }
+
+ private:
+ /// Matrix, kept as a member only for solving successive linear equation
+ /// that preconditioner matrix may vary.
+ PETScMatrix &_A;
+
+ KSP _solver; ///< Solver type.
+ PC _pc; ///< Preconditioner type.
+};
+
+} // end namespace
+#endif
+
diff --git a/MathLib/LinAlg/PETSc/PETScMatrix.cpp b/MathLib/LinAlg/PETSc/PETScMatrix.cpp
index 2cafefa044544c331b6332991d91eae0c73d4cb9..8e706128cb03fc60d6e79f3a2168532c27c4b29e 100644
--- a/MathLib/LinAlg/PETSc/PETScMatrix.cpp
+++ b/MathLib/LinAlg/PETSc/PETScMatrix.cpp
@@ -75,7 +75,7 @@ void PETScMatrix::viewer(const std::string &file_name, const PetscViewerFormat v
MatView(_A,viewer);
// This preprocessor is only for debugging, e.g. dump the matrix and exit the program.
-//#define EXIT_TEST
+//#define EXIT_TEST
#ifdef EXIT_TEST
MatDestroy(&_A);
PetscFinalize();
@@ -91,7 +91,7 @@ void PETScMatrix::create(const PetscInt d_nz, const PetscInt o_nz)
MatSetFromOptions(_A);
- // for a dense matrix: MatSeqAIJSetPreallocation(_A, d_nz, PETSC_NULL);
+ MatSeqAIJSetPreallocation(_A, d_nz, PETSC_NULL);
MatMPIAIJSetPreallocation(_A, d_nz, PETSC_NULL, o_nz, PETSC_NULL);
MatGetOwnershipRange(_A, &_start_rank, &_end_rank);
diff --git a/MathLib/LinAlg/PETSc/PETScTools.cpp b/MathLib/LinAlg/PETSc/PETScTools.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..aea99c360604df6bd2d443733072e37f69b99010
--- /dev/null
+++ b/MathLib/LinAlg/PETSc/PETScTools.cpp
@@ -0,0 +1,41 @@
+/*!
+ \file PETScTools.cpp
+ \brief Definition of a function related to PETSc solver interface to assign
+ the Dirichlet boundary conditions.
+
+ \author Wenqing Wang
+ \version
+ \date Nov 2011 - Sep 2013
+
+ \copyright
+ Copyright (c) 2013, OpenGeoSys Community (http://www.opengeosys.org)
+ Distributed under a Modified BSD License.
+ See accompanying file LICENSE.txt or
+ http://www.opengeosys.org/project/license
+*/
+
+#include "PETScTools.h"
+
+namespace MathLib
+{
+
+void applyKnownSolution(PETScMatrix &A, PETScVector &b, PETScVector &x,
+ const std::vector<PetscInt> &vec_knownX_id,
+ const std::vector<PetscScalar> &vec_knownX_x)
+{
+ A.setRowsColumnsZero(vec_knownX_id);
+ A.finalizeAssembly();
+
+ if(vec_knownX_id.size() > 0)
+ {
+ x.set(vec_knownX_id, vec_knownX_x);
+ b.set(vec_knownX_id, vec_knownX_x);
+ }
+
+ x.finalizeAssembly();
+ b.finalizeAssembly();
+}
+
+} // end of namespace MathLib
+
+
diff --git a/MathLib/LinAlg/PETSc/PETScTools.h b/MathLib/LinAlg/PETSc/PETScTools.h
new file mode 100644
index 0000000000000000000000000000000000000000..aa9011719c2e01918368a4bfee54a7feb526ce3d
--- /dev/null
+++ b/MathLib/LinAlg/PETSc/PETScTools.h
@@ -0,0 +1,43 @@
+/*!
+ \file PETScTools.h
+ \brief Declaration of a function related to PETSc solver interface to assign
+ the Dirichlet boundary conditions.
+
+ \author Wenqing Wang
+ \version
+ \date Nov 2011 - Sep 2013
+
+ \copyright
+ Copyright (c) 2013, OpenGeoSys Community (http://www.opengeosys.org)
+ Distributed under a Modified BSD License.
+ See accompanying file LICENSE.txt or
+ http://www.opengeosys.org/project/license
+*/
+
+#ifndef PETSCTOOLS_H_
+#define PETSCTOOLS_H_
+
+#include <vector>
+
+#include "PETScMatrix.h"
+#include "PETScVector.h"
+
+namespace MathLib
+{
+/*!
+ \brief apply known solutions to a system of linear equations
+
+ \param A Coefficient matrix
+ \param b RHS vector
+ \param vec_knownX_id a vector of known solution entry IDs
+ \param vec_knownX_x a vector of known solutions
+ \param penalty_scaling value for scaling some matrix and right hand side
+ entries to enforce some conditions
+*/
+void applyKnownSolution(PETScMatrix &A, PETScVector &b, PETScVector &x,
+ const std::vector<PetscInt> &vec_knownX_id,
+ const std::vector<PetscScalar> &vec_knownX_x);
+} // end of namespace MathLib
+
+#endif //end of PETSCTOOLS_H_
+
diff --git a/MathLib/LinAlg/Solvers/GaussAlgorithm-impl.h b/MathLib/LinAlg/Solvers/GaussAlgorithm-impl.h
index 22a31166ca7428bf1a5d83071d93cfc8a8812c68..25a6fa485ff70ecda62e9d6b036e46b6493df80b 100644
--- a/MathLib/LinAlg/Solvers/GaussAlgorithm-impl.h
+++ b/MathLib/LinAlg/Solvers/GaussAlgorithm-impl.h
@@ -32,7 +32,7 @@ GaussAlgorithm<MAT_T, VEC_T>::GaussAlgorithm(MAT_T &A,
_perm[k] = k;
for (i=k+1; i<nr; i++) {
if (std::abs(_mat(i,k)) > t) {
- t = _mat(i,k);
+ t = std::abs(_mat(i,k));
_perm[k] = i;
}
}
diff --git a/MathLib/MathTools.cpp b/MathLib/MathTools.cpp
index 651fa75369b10a4218753378150339db1694d346..a086660e996d0f8a01256e1c90eb0ba8a49489c7 100644
--- a/MathLib/MathTools.cpp
+++ b/MathLib/MathTools.cpp
@@ -62,23 +62,6 @@ double getAngle (const double p0[3], const double p1[3], const double p2[3])
return acos (scalarProduct<double,3> (v0,v1) / (sqrt(scalarProduct<double,3>(v0,v0)) * sqrt(scalarProduct<double,3>(v1,v1))));
}
-double calcTriangleArea(const double* p0, const double* p1, const double* p2)
-{
- const double u0 (p2[0] - p0[0]);
- const double u1 (p2[1] - p0[1]);
- const double u2 (p2[2] - p0[2]);
-
- const double v0 (p1[0] - p0[0]);
- const double v1 (p1[1] - p0[1]);
- const double v2 (p1[2] - p0[2]);
-
- const double z0 (u1 * v2 - u2 * v1);
- const double z1 (u2 * v0 - u0 * v2);
- const double z2 (u0 * v1 - u1 * v0);
-
- return 0.5 * sqrt(z0 * z0 + z1 * z1 + z2 * z2);
-}
-
double calcTetrahedronVolume(const double* x1, const double* x2, const double* x3, const double* x4)
{
return fabs((x1[0] - x4[0]) * ((x2[1] - x4[1]) * (x3[2] - x4[2]) - (x2[2] - x4[2]) * (x3[1] - x4[1]))
diff --git a/MathLib/MathTools.h b/MathLib/MathTools.h
index 585d141b14c405b6143eeb70d2ac50a59fc559ab..8c78efca61690d9c1f4cc4dca30118855648d0af 100644
--- a/MathLib/MathTools.h
+++ b/MathLib/MathTools.h
@@ -148,12 +148,6 @@ float normalize(float min, float max, float val);
*/
double getAngle (const double p0[3], const double p1[3], const double p2[3]);
-/**
- * Calculates the area of a triangle.
- * The formula is A=.5*|u x v|, i.e. half of the area of the parallelogram specified by u=p0->p1 and v=p0->p2.
- */
-double calcTriangleArea(const double* p0, const double* p1, const double* p2);
-
/**
* Calculates the volume of a tetrahedron.
* The formula is V=1/6*|a(b x c)| with a=x1->x2, b=x1->x3 and c=x1->x4.
diff --git a/MeshGeoToolsLib/AppendLinesAlongPolyline.cpp b/MeshGeoToolsLib/AppendLinesAlongPolyline.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..074983358b340bb288879f4aaa6e12df62bc1e4a
--- /dev/null
+++ b/MeshGeoToolsLib/AppendLinesAlongPolyline.cpp
@@ -0,0 +1,69 @@
+
+/**
+ * @copyright
+ * Copyright (c) 2013, OpenGeoSys Community (http://www.opengeosys.org)
+ * Distributed under a Modified BSD License.
+ * See accompanying file LICENSE.txt or
+ * http://www.opengeosys.org/LICENSE.txt
+ */
+
+#include "AppendLinesAlongPolyline.h"
+
+// ThirdParty/logog
+#include "logog/include/logog.hpp"
+
+// MeshLib
+#include "Mesh.h"
+#include "Node.h"
+#include "Elements/Line.h"
+#include "Elements/Element.h"
+#include "MeshEnums.h"
+#include "MeshEditing/DuplicateMeshComponents.h"
+
+#include "MeshGeoToolsLib/MeshNodesAlongPolyline.h"
+
+namespace MeshGeoToolsLib
+{
+
+MeshLib::Mesh* appendLinesAlongPolylines(const MeshLib::Mesh &mesh, const GeoLib::PolylineVec &ply_vec)
+{
+
+ // copy existing nodes and elements
+ std::vector<MeshLib::Node*> vec_new_nodes = MeshLib::copyNodeVector(mesh.getNodes());
+ std::vector<MeshLib::Element*> vec_new_eles = MeshLib::copyElementVector(mesh.getElements(), vec_new_nodes);
+ unsigned max_matID = 0;
+ for (auto e : vec_new_eles)
+ max_matID = std::max(max_matID, e->getValue());
+
+ const size_t n_ply (ply_vec.size());
+ // for each polyline
+ for (size_t k(0); k < n_ply; k++)
+ {
+ const GeoLib::Polyline* ply = (*ply_vec.getVector())[k];
+
+ // search nodes on the polyline
+ MeshGeoToolsLib::MeshNodesAlongPolyline mshNodesAlongPoly(mesh, *ply, mesh.getMinEdgeLength()*0.5);
+ auto &vec_nodes_on_ply = mshNodesAlongPoly.getNodeIDs();
+ if (vec_nodes_on_ply.empty()) {
+ std::string ply_name;
+ ply_vec.getNameOfElementByID(k, ply_name);
+ INFO("No nodes found on polyline %s", ply_name.c_str());
+ continue;
+ }
+
+ // add line elements
+ for (size_t i=0; i<vec_nodes_on_ply.size()-1; i++) {
+ std::array<MeshLib::Node*, 2> element_nodes;
+ element_nodes[0] = vec_new_nodes[vec_nodes_on_ply[i]];
+ element_nodes[1] = vec_new_nodes[vec_nodes_on_ply[i+1]];
+ vec_new_eles.push_back(new MeshLib::Line(element_nodes, max_matID+k+1));
+ }
+ }
+
+ // generate a mesh
+ const std::string new_mesh_name = mesh.getName() + "_with_lines";
+ return new MeshLib::Mesh(new_mesh_name, vec_new_nodes, vec_new_eles);
+}
+
+} // MeshGeoToolsLib
+
diff --git a/MeshGeoToolsLib/AppendLinesAlongPolyline.h b/MeshGeoToolsLib/AppendLinesAlongPolyline.h
new file mode 100644
index 0000000000000000000000000000000000000000..522d05914e004549375cd72fcd1255d95d1f79cc
--- /dev/null
+++ b/MeshGeoToolsLib/AppendLinesAlongPolyline.h
@@ -0,0 +1,40 @@
+/**
+ * @copyright
+ * Copyright (c) 2013, OpenGeoSys Community (http://www.opengeosys.org)
+ * Distributed under a Modified BSD License.
+ * See accompanying file LICENSE.txt or
+ * http://www.opengeosys.org/LICENSE.txt
+ */
+
+#ifndef APPENDLINESALONGPOLYLINES_H_
+#define APPENDLINESALONGPOLYLINES_H_
+
+// GeoLib
+#include "PolylineVec.h"
+
+namespace MeshLib
+{
+class Mesh;
+}
+
+namespace MeshGeoToolsLib
+{
+
+/**
+ * Add line elements to a copy of a given mesh
+ *
+ * The function creates line elements from nodes located along user-provided polylines.
+ * New elements will have a distinct material ID for each polyline.
+ *
+ * \remark The function allows creation of duplicated line elements.
+ * \remark Line elements may not be placed along edges of existing elements.
+ *
+ * @param mesh original mesh
+ * @param ply_vec polyline vector whose nodes are used to create line elements
+ * @return a new mesh which is copied from a given mesh and additionally includes line elements
+ */
+MeshLib::Mesh* appendLinesAlongPolylines(const MeshLib::Mesh &mesh, const GeoLib::PolylineVec &ply_vec);
+
+}
+
+#endif
diff --git a/MeshLib/Elements/TemplateQuad-impl.h b/MeshLib/Elements/TemplateQuad-impl.h
index afd624856b42bf684cc3ca4911314694970c72c8..fff572636170247c3932311c97c6387247f765fb 100644
--- a/MeshLib/Elements/TemplateQuad-impl.h
+++ b/MeshLib/Elements/TemplateQuad-impl.h
@@ -90,8 +90,8 @@ TemplateQuad<NNODES,CELLQUADTYPE>::~TemplateQuad()
template <unsigned NNODES, CellType CELLQUADTYPE>
double TemplateQuad<NNODES,CELLQUADTYPE>::computeVolume()
{
- return MathLib::calcTriangleArea(_nodes[0]->getCoords(), _nodes[1]->getCoords(), _nodes[2]->getCoords())
- + MathLib::calcTriangleArea(_nodes[2]->getCoords(), _nodes[3]->getCoords(), _nodes[0]->getCoords());
+ return GeoLib::calcTriangleArea(*_nodes[0], *_nodes[1], *_nodes[2])
+ + GeoLib::calcTriangleArea(*_nodes[2], *_nodes[3], *_nodes[0]);
}
template <unsigned NNODES, CellType CELLQUADTYPE>
diff --git a/MeshLib/Elements/TemplateTri.h b/MeshLib/Elements/TemplateTri.h
index a0c05f27d38dd4ccb0cf4a630e3eb3d951b56cb4..c91290193e21420ee7fa232d2bc5aab822f78fa7 100644
--- a/MeshLib/Elements/TemplateTri.h
+++ b/MeshLib/Elements/TemplateTri.h
@@ -22,6 +22,7 @@
#include "MeshEnums.h"
#include "MathTools.h"
+#include "AnalyticalGeometry.h"
namespace MeshLib {
@@ -121,7 +122,7 @@ protected:
/// Calculates the area of the triangle by returning half of the area of the corresponding parallelogram.
double computeVolume()
{
- return MathLib::calcTriangleArea(_nodes[0]->getCoords(), _nodes[1]->getCoords(), _nodes[2]->getCoords());
+ return GeoLib::calcTriangleArea(*_nodes[0], *_nodes[1], *_nodes[2]);
}
protected:
diff --git a/MeshLib/MeshGenerators/MeshGenerator.cpp b/MeshLib/MeshGenerators/MeshGenerator.cpp
index cd6576856bab399266c7d8e1183e391d04aea06c..2f30171f6780cd9534df54e57aafd33bdef24598 100644
--- a/MeshLib/MeshGenerators/MeshGenerator.cpp
+++ b/MeshLib/MeshGenerators/MeshGenerator.cpp
@@ -18,6 +18,7 @@
#include "Elements/Line.h"
#include "Elements/Quad.h"
#include "Elements/Hex.h"
+#include "Elements/Tri.h"
#include <vector>
@@ -174,4 +175,57 @@ Mesh* MeshGenerator::generateRegularHexMesh(const unsigned n_x_cells,
return new Mesh(mesh_name, nodes, elements);
}
+Mesh* MeshGenerator::generateRegularTriMesh(
+ const double length,
+ const std::size_t subdivision,
+ const GeoLib::Point& origin)
+{
+ return generateRegularTriMesh(subdivision, subdivision, length/subdivision, origin);
+}
+
+Mesh* MeshGenerator::generateRegularTriMesh(const unsigned n_x_cells,
+ const unsigned n_y_cells,
+ const double cell_size,
+ GeoLib::Point const& origin,
+ std::string const& mesh_name)
+{
+ //nodes
+ const unsigned n_x_nodes (n_x_cells+1);
+ const unsigned n_y_nodes (n_y_cells+1);
+ std::vector<Node*> nodes;
+ nodes.reserve(n_x_nodes * n_y_nodes);
+
+ for (std::size_t i = 0; i < n_y_nodes; i++)
+ {
+ const double y_offset (origin[1] + cell_size * i);
+ for (std::size_t j = 0; j < n_x_nodes; j++)
+ nodes.push_back (new Node(origin[0] + cell_size * j, y_offset, origin[2]));
+ }
+
+ //elements
+ std::vector<Element*> elements;
+ elements.reserve(n_x_cells * n_y_cells * 2);
+
+ for (std::size_t j = 0; j < n_y_cells; j++)
+ {
+ const std::size_t offset_y1 = j * n_x_nodes;
+ const std::size_t offset_y2 = (j + 1) * n_x_nodes;
+ for (std::size_t k = 0; k < n_x_cells; k++)
+ {
+ std::array<Node*, 3> element1_nodes;
+ element1_nodes[0] = nodes[offset_y1 + k];
+ element1_nodes[1] = nodes[offset_y2 + k + 1];
+ element1_nodes[2] = nodes[offset_y2 + k];
+ elements.push_back (new Tri(element1_nodes));
+ std::array<Node*, 3> element2_nodes;
+ element2_nodes[0] = nodes[offset_y1 + k];
+ element2_nodes[1] = nodes[offset_y1 + k + 1];
+ element2_nodes[2] = nodes[offset_y2 + k + 1];
+ elements.push_back (new Tri(element2_nodes));
+ }
+ }
+
+ return new Mesh(mesh_name, nodes, elements);
+}
+
}
diff --git a/MeshLib/MeshGenerators/MeshGenerator.h b/MeshLib/MeshGenerators/MeshGenerator.h
index d9f72dc833e368023cdbbe5042c76fede33be0fb..11946d2711e80af1cc6a385e23cd573b43b384f2 100644
--- a/MeshLib/MeshGenerators/MeshGenerator.h
+++ b/MeshLib/MeshGenerators/MeshGenerator.h
@@ -103,6 +103,34 @@ Mesh* generateRegularHexMesh(const unsigned n_x_cells,
GeoLib::Point const& origin = GeoLib::ORIGIN,
std::string const& mesh_name = "mesh");
+/**
+ * Generate a regular 2D Triangle-Element mesh. The mesh is generated in the
+ * x-y-plane.
+ *
+ * \param length Mesh's dimensions in x- and y-directions.
+ * \param subdivision Number of subdivisions.
+ * \param origin Optional mesh's origin (the lower left corner) with GeoLib::ORIGIN default.
+ */
+Mesh* generateRegularTriMesh(const double length,
+ const std::size_t subdivision,
+ GeoLib::Point const& origin = GeoLib::ORIGIN);
+
+/**
+ * Generate a regular 2D Triangle-Element mesh. The mesh is generated in the
+ * x-y-plane.
+ *
+ * \param n_x_cells Number of cells in x-direction.
+ * \param n_y_cells Number of cells in y-direction.
+ * \param cell_size Edge length of two equal sides of isosceles triangles
+ * \param origin Optional mesh's origin (the lower left corner) with GeoLib::ORIGIN default.
+ * \param mesh_name Name of the new mesh.
+ */
+Mesh* generateRegularTriMesh(const unsigned n_x_cells,
+ const unsigned n_y_cells,
+ const double cell_size,
+ GeoLib::Point const& origin = GeoLib::ORIGIN,
+ std::string const& mesh_name = "mesh");
+
} //MeshGenerator
} //MeshLib
diff --git a/NumLib/Function/LinearInterpolationOnSurface.cpp b/NumLib/Function/LinearInterpolationOnSurface.cpp
index d71716a7b856ac5400e1ea1ccacd5f5b105e328d..18aa67f09ae1f3a19cf025b85722bd54756d1429 100644
--- a/NumLib/Function/LinearInterpolationOnSurface.cpp
+++ b/NumLib/Function/LinearInterpolationOnSurface.cpp
@@ -70,10 +70,10 @@ double LinearInterpolationOnSurface::interpolateInTri(const GeoLib::Triangle &tr
std::vector<GeoLib::Point*> p_pnts = {{&pnts[0], &pnts[1], &pnts[2]}};
GeoLib::rotatePointsToXY(p_pnts);
- const double* v1 = pnts[0].getCoords();
- const double* v2 = pnts[1].getCoords();
- const double* v3 = pnts[2].getCoords();
- const double area = MathLib::calcTriangleArea(v1, v2, v3);
+ GeoLib::Point const& v1(pnts[0]);
+ GeoLib::Point const& v2(pnts[1]);
+ GeoLib::Point const& v3(pnts[2]);
+ const double area = GeoLib::calcTriangleArea(v1, v2, v3);
if (area==.0) {
// take average if all points have the same coordinates
diff --git a/SimpleTests/MatrixTests/CMakeLists.txt b/SimpleTests/MatrixTests/CMakeLists.txt
index 9a67ca86c403c8102bd999301f383b804a1b0508..00ae009881fe63fed08f601cf7336aab346607fa 100644
--- a/SimpleTests/MatrixTests/CMakeLists.txt
+++ b/SimpleTests/MatrixTests/CMakeLists.txt
@@ -101,3 +101,15 @@ TARGET_LINK_LIBRARIES ( MatTestRemoveRowsCols
MathLib
)
+ADD_EXECUTABLE( DenseGaussEliminationChecker
+ DenseGaussEliminationChecker.cpp
+ ${SOURCES}
+ ${HEADERS}
+)
+SET_TARGET_PROPERTIES(DenseGaussEliminationChecker PROPERTIES FOLDER SimpleTests)
+TARGET_LINK_LIBRARIES(DenseGaussEliminationChecker
+ logog
+ BaseLib
+ MathLib
+)
+
diff --git a/SimpleTests/MatrixTests/DenseGaussEliminationChecker.cpp b/SimpleTests/MatrixTests/DenseGaussEliminationChecker.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..d8b5fccf9c9e8387a812e767dfcc1f44dff20525
--- /dev/null
+++ b/SimpleTests/MatrixTests/DenseGaussEliminationChecker.cpp
@@ -0,0 +1,88 @@
+/**
+ * \date 2014-06-11
+ * \brief Implementation of tests.
+ *
+ * \copyright
+ * Copyright (c) 2014, OpenGeoSys Community (http://www.opengeosys.org)
+ * Distributed under a Modified BSD License.
+ * See accompanying file LICENSE.txt or
+ * http://www.opengeosys.org/project/license
+ *
+ */
+
+#include <fstream>
+#include <sstream>
+
+// BaseLib/tclap
+#include "tclap/CmdLine.h"
+// ThirdParty/logog
+#include "logog/include/logog.hpp"
+#include "logog/include/formatter.hpp"
+
+// BaseLib
+#include "LogogSimpleFormatter.h"
+
+// MathLib
+#include "LinAlg/Dense/DenseMatrix.h"
+#include "LinAlg/Solvers/GaussAlgorithm.h"
+
+int main(int argc, char *argv[])
+{
+ LOGOG_INITIALIZE();
+ BaseLib::LogogSimpleFormatter *custom_format (new BaseLib::LogogSimpleFormatter);
+ logog::Cout *logogCout(new logog::Cout);
+ logogCout->SetFormatter(*custom_format);
+
+ TCLAP::CmdLine cmd("Simple direct matrix solver test.\n\
+ It consists of the following steps:\n\
+ (1) Read a matrix A from ascii format\n\
+ (2) Set all entries of a vector x to one and compute b = A * x\n\
+ (3) Solve the system of linear equations -> result have to be (1,...,1)", ' ', "0.1");
+ TCLAP::ValueArg<std::string> matrix_arg("m", "matrix", "input matrix file (ascii format)", true, "", "string");
+ cmd.add( matrix_arg );
+ cmd.parse( argc, argv );
+
+ // *** reading dense matrix in ascii format from file
+ std::string const fname_mat(matrix_arg.getValue());
+ std::ifstream in(fname_mat.c_str());
+ if (!in) {
+ ERR("error reading matrix from %s", fname_mat.c_str());
+ return -1;
+ }
+ INFO("reading matrix from %s ...", fname_mat.c_str());
+ std::size_t n_rows(0), n_cols(0);
+ in >> n_rows;
+ in >> n_cols;
+ MathLib::DenseMatrix<double, std::size_t> mat(n_rows, n_cols);
+ for (std::size_t i(0); i<mat.getNRows(); ++i) {
+ for (std::size_t j(0); j<mat.getNCols(); ++j) {
+ in >> mat(i,j);
+ }
+ }
+ {
+ std::stringstream stream;
+ stream << mat;
+ INFO("read matrix:\n%s", stream.str().c_str());
+ }
+
+ std::vector<double> x(n_cols,1.0), b;
+ b.resize(n_rows);
+ b = mat * x;
+
+ MathLib::GaussAlgorithm<MathLib::DenseMatrix<double, std::size_t>> gauss(mat);
+ gauss.solve(b);
+
+ {
+ std::stringstream stream;
+ std::copy(b.begin(), b.end(), std::ostream_iterator<double>(stream, " "));
+ stream << std::endl;
+ INFO("solution vector:\n%s", stream.str().c_str());
+ }
+
+ delete custom_format;
+ delete logogCout;
+ LOGOG_SHUTDOWN();
+
+ return 0;
+}
+
diff --git a/Tests/CMakeLists.txt b/Tests/CMakeLists.txt
index ab1362afa47e378fe2adeeb3c8149514d3f70899..3e5cc42b013d09be328a19dcdd635192f21011f3 100644
--- a/Tests/CMakeLists.txt
+++ b/Tests/CMakeLists.txt
@@ -78,7 +78,11 @@ ENDIF()
# This should override CTest's predefined test-target but it does not
IF (OGS_USE_PETSC)
ADD_CUSTOM_TARGET(test
- mpirun -np 3 $<TARGET_FILE:testrunner>
+ mpirun -np 1 $<TARGET_FILE:testrunner> --gtest_filter=-MPITest*
+ DEPENDS testrunner
+ )
+ ADD_CUSTOM_TARGET(test_mpi
+ mpirun -np 3 $<TARGET_FILE:testrunner> --gtest_filter=MPITest*
DEPENDS testrunner
)
ELSE ()
diff --git a/Tests/GeoLib/TestComputeAndInsertAllIntersectionPoints.cpp b/Tests/GeoLib/TestComputeAndInsertAllIntersectionPoints.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..1495ed8a33c0ed33fc6c3993566ae64df37c962f
--- /dev/null
+++ b/Tests/GeoLib/TestComputeAndInsertAllIntersectionPoints.cpp
@@ -0,0 +1,108 @@
+/**
+ * @file TestComputeAndInsertAllIntersectionPoints.cpp
+ * @date 2014-05-09
+ *
+ * @copyright
+ * Copyright (c) 2014, OpenGeoSys Community (http://www.opengeosys.org)
+ * Distributed under a Modified BSD License.
+ * See accompanying file LICENSE.txt or
+ * http://www.opengeosys.org/LICENSE.txt
+ */
+
+#include <ctime>
+#include <tuple>
+
+#include "gtest/gtest.h"
+
+#include "Point.h"
+#include "Polyline.h"
+#include "GEOObjects.h"
+#include "AnalyticalGeometry.h"
+
+TEST(GeoLib, TestComputeAndInsertAllIntersectionPoints)
+{
+ // *** insert points in vector
+ std::vector<GeoLib::Point*> *pnts(new std::vector<GeoLib::Point*>);
+ pnts->push_back(new GeoLib::Point(0.0,0.0,0.0));
+ pnts->push_back(new GeoLib::Point(11.0,0.0,0.0));
+
+ pnts->push_back(new GeoLib::Point(0.0,1.0,0.0));
+ pnts->push_back(new GeoLib::Point(1.0,-1.0,0.0));
+ pnts->push_back(new GeoLib::Point(2.0, 1.0,0.0));
+ pnts->push_back(new GeoLib::Point(3.0,-1.0,0.0));
+ pnts->push_back(new GeoLib::Point(4.0, 1.0,0.0));
+ pnts->push_back(new GeoLib::Point(5.0,-1.0,0.0));
+ pnts->push_back(new GeoLib::Point(6.0, 1.0,0.0));
+ pnts->push_back(new GeoLib::Point(7.0,-1.0,0.0));
+ pnts->push_back(new GeoLib::Point(8.0, 1.0,0.0));
+ pnts->push_back(new GeoLib::Point(9.0,-1.0,0.0));
+ pnts->push_back(new GeoLib::Point(10.0, 1.0,0.0));
+ pnts->push_back(new GeoLib::Point(11.0,-1.0,0.0));
+
+ GeoLib::GEOObjects geo_objs;
+ std::string geo_name("TestGeometry");
+ geo_objs.addPointVec(pnts, geo_name);
+
+ // *** create polylines
+ GeoLib::Polyline* ply0(new GeoLib::Polyline(*pnts));
+ ply0->addPoint(0);
+ ply0->addPoint(1);
+ GeoLib::Polyline* ply1(new GeoLib::Polyline(*pnts));
+ for (std::size_t k(2); k<14; k++)
+ ply1->addPoint(k);
+ std::vector<GeoLib::Polyline*>* plys(new std::vector<GeoLib::Polyline*>);
+ plys->push_back(ply0);
+ plys->push_back(ply1);
+
+ GeoLib::PointVec &pnt_vec(*(const_cast<GeoLib::PointVec*>(geo_objs.getPointVecObj(geo_name))));
+ GeoLib::computeAndInsertAllIntersectionPoints(pnt_vec, *plys);
+
+ ASSERT_EQ(25u, pnt_vec.size());
+ ASSERT_EQ(13u, ply0->getNumberOfPoints());
+ ASSERT_EQ(23u, ply1->getNumberOfPoints());
+
+ // check correct order of points in ply0
+ EXPECT_EQ(0u, ply0->getPointID(0));
+ EXPECT_EQ(14u, ply0->getPointID(1));
+ EXPECT_EQ(15u, ply0->getPointID(2));
+ EXPECT_EQ(16u, ply0->getPointID(3));
+ EXPECT_EQ(17u, ply0->getPointID(4));
+ EXPECT_EQ(18u, ply0->getPointID(5));
+ EXPECT_EQ(19u, ply0->getPointID(6));
+ EXPECT_EQ(20u, ply0->getPointID(7));
+ EXPECT_EQ(21u, ply0->getPointID(8));
+ EXPECT_EQ(22u, ply0->getPointID(9));
+ EXPECT_EQ(23u, ply0->getPointID(10));
+ EXPECT_EQ(24u, ply0->getPointID(11));
+ EXPECT_EQ(1u, ply0->getPointID(12));
+
+ // check correct order of points in ply1
+ EXPECT_EQ(2u, ply1->getPointID(0));
+ EXPECT_EQ(14u, ply1->getPointID(1));
+ EXPECT_EQ(3u, ply1->getPointID(2));
+ EXPECT_EQ(15u, ply1->getPointID(3));
+ EXPECT_EQ(4u, ply1->getPointID(4));
+ EXPECT_EQ(16u, ply1->getPointID(5));
+ EXPECT_EQ(5u, ply1->getPointID(6));
+ EXPECT_EQ(17u, ply1->getPointID(7));
+ EXPECT_EQ(6u, ply1->getPointID(8));
+ EXPECT_EQ(18u, ply1->getPointID(9));
+ EXPECT_EQ(7u, ply1->getPointID(10));
+ EXPECT_EQ(19u, ply1->getPointID(11));
+ EXPECT_EQ(8u, ply1->getPointID(12));
+ EXPECT_EQ(20u, ply1->getPointID(13));
+ EXPECT_EQ(9u, ply1->getPointID(14));
+ EXPECT_EQ(21u, ply1->getPointID(15));
+ EXPECT_EQ(10u, ply1->getPointID(16));
+ EXPECT_EQ(22u, ply1->getPointID(17));
+ EXPECT_EQ(11u, ply1->getPointID(18));
+ EXPECT_EQ(23u, ply1->getPointID(19));
+ EXPECT_EQ(12u, ply1->getPointID(20));
+ EXPECT_EQ(24u, ply1->getPointID(21));
+ EXPECT_EQ(13u, ply1->getPointID(22));
+
+ delete plys;
+ delete ply1;
+ delete ply0;
+}
+
diff --git a/Tests/GeoLib/TestIsPointInTriangle.cpp b/Tests/GeoLib/TestIsPointInTriangle.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..ce90557fa3a5b82c894d91f512649f7ed03dc1cf
--- /dev/null
+++ b/Tests/GeoLib/TestIsPointInTriangle.cpp
@@ -0,0 +1,60 @@
+/**
+ * @file TestIsPointInTriangle.cpp
+ * @date 2014-05-27
+ *
+ * @copyright
+ * Copyright (c) 2014, OpenGeoSys Community (http://www.opengeosys.org)
+ * Distributed under a Modified BSD License.
+ * See accompanying file LICENSE.txt or
+ * http://www.opengeosys.org/LICENSE.txt
+ */
+
+#include "gtest/gtest.h"
+
+#include "AnalyticalGeometry.h"
+
+TEST(GeoLib, IsPointInTriangle)
+{
+ GeoLib::Point const a(0.0, 0.0, 0.0);
+ GeoLib::Point const b(100.0, 0.0, 0.0);
+ GeoLib::Point const c(0.0, 100.0, 0.0);
+
+ // check point on corner points of triangle
+ GeoLib::Point q(a);
+ EXPECT_TRUE(GeoLib::isPointInTriangle(q, a, b, c));
+ q = b;
+ EXPECT_TRUE(GeoLib::isPointInTriangle(q, a, b, c));
+ q = c;
+ EXPECT_TRUE(GeoLib::isPointInTriangle(q, a, b, c));
+
+ // check points on edges of triangle
+ q = GeoLib::Point(0.5*(b[0]+a[0]), 0.5*(b[1]+a[1]), 0.5*(b[2]+a[2]));
+ EXPECT_TRUE(GeoLib::isPointInTriangle(q, a, b, c));
+ q = GeoLib::Point(0.5*(c[0]+a[0]), 0.5*(c[1]+a[1]), 0.5*(c[2]+a[2]));
+ EXPECT_TRUE(GeoLib::isPointInTriangle(q, a, b, c));
+ q = GeoLib::Point(0.5*(c[0]+b[0]), 0.5*(c[1]+b[1]), 0.5*(c[2]+b[2]));
+ EXPECT_TRUE(GeoLib::isPointInTriangle(q, a, b, c));
+
+ // check points inside
+ q = GeoLib::Point (0.1, 0.1, 0.0);
+ EXPECT_TRUE(GeoLib::isPointInTriangle(q, a, b, c));
+ q = GeoLib::Point (0.1, 0.1, 1e-10);
+ EXPECT_TRUE(GeoLib::isPointInTriangle(q, a, b, c));
+
+ // check points outside
+ q = GeoLib::Point (-0.1, 0.1, 0.0);
+ EXPECT_FALSE(GeoLib::isPointInTriangle(q, a, b, c));
+ q = GeoLib::Point (0.1, 0.1, 0.0005);
+ EXPECT_FALSE(GeoLib::isPointInTriangle(q, a, b, c)); // default eps from float
+ q = GeoLib::Point (0.1, 0.1, 0.0001);
+ EXPECT_FALSE(GeoLib::isPointInTriangle(q, a, b, c, std::numeric_limits<double>::epsilon()));
+ q = GeoLib::Point (0.1, 0.1, 0.000001);
+ EXPECT_FALSE(GeoLib::isPointInTriangle(q, a, b, c, std::numeric_limits<double>::epsilon()));
+ q = GeoLib::Point (0.1, 0.1, 1e-7);
+ EXPECT_FALSE(GeoLib::isPointInTriangle(q, a, b, c, std::numeric_limits<double>::epsilon()));
+ q = GeoLib::Point (0.1, 0.1, 0.001);
+ EXPECT_FALSE(GeoLib::isPointInTriangle(q, a, b, c));
+ q = GeoLib::Point (0.1, 0.1, 0.1);
+ EXPECT_FALSE(GeoLib::isPointInTriangle(q, a, b, c));
+}
+
diff --git a/Tests/GeoLib/TestPolygon.cpp b/Tests/GeoLib/TestPolygon.cpp
index 47444e0ca5fc63445121fe6d00ea8c9227ef04bf..eeaa7b8155c8854a2c57924b481e705f3c870b2f 100644
--- a/Tests/GeoLib/TestPolygon.cpp
+++ b/Tests/GeoLib/TestPolygon.cpp
@@ -20,32 +20,34 @@ using namespace GeoLib;
/**
* Polygon:
- * 3
- * / \
- * / \
- * 2 4
- * | |
- * | |
- * 1 5
- * \ /
- * \ /
- * 0
+ * 2 4 6
+ * |\ / \ /|
+ * | \ / \ / |
+ * 1 3 5 7
+ * \ /
+ * \ /
+ * \ /
+ * \ /
+ * \ /
+ * 0
*/
-class IsPntInPolygonTest : public testing::Test
+class PolygonTest : public testing::Test
{
public:
- IsPntInPolygonTest() :
+ PolygonTest() :
_polygon(nullptr)
{
// create points and construct polygon
- _pnts.push_back(new Point(1.0,0.0,0.0));
- _pnts.push_back(new Point(0.0,1.0,0.0));
- _pnts.push_back(new Point(0.0,2.0,0.0));
- _pnts.push_back(new Point(1.0,3.0,0.0));
- _pnts.push_back(new Point(2.0,2.0,0.0));
- _pnts.push_back(new Point(2.0,1.0,0.0));
+ _pnts.push_back(new Point( 0.0, 0.0,0.0)); // 0
+ _pnts.push_back(new Point(-2.0, 2.0,0.0)); // 1
+ _pnts.push_back(new Point(-2.0, 4.0,0.0)); // 2
+ _pnts.push_back(new Point(-1.0, 2.0,0.0)); // 3
+ _pnts.push_back(new Point( 0.0, 4.0,0.0)); // 4
+ _pnts.push_back(new Point( 1.0, 2.0,0.0)); // 5
+ _pnts.push_back(new Point( 2.0, 4.0,0.0)); // 6
+ _pnts.push_back(new Point( 2.0, 2.0,0.0)); // 7
// create closed polyline
Polyline ply(_pnts);
@@ -55,13 +57,15 @@ public:
ply.addPoint(3);
ply.addPoint(4);
ply.addPoint(5);
+ ply.addPoint(6);
+ ply.addPoint(7);
ply.addPoint(0);
// create polygon
_polygon = new Polygon(ply);
}
- ~IsPntInPolygonTest()
+ ~PolygonTest()
{
delete _polygon;
for (std::size_t k(0); k<_pnts.size(); k++)
@@ -73,13 +77,13 @@ protected:
Polygon *_polygon;
};
-TEST_F(IsPntInPolygonTest, CheckCorners)
+TEST_F(PolygonTest, isPntInPolygonCheckCorners)
{
for (std::size_t k(0); k<_pnts.size(); k++)
EXPECT_TRUE(_polygon->isPntInPolygon(*_pnts[k]));
}
-TEST_F(IsPntInPolygonTest, CheckPointsRestOnPolygonEdges)
+TEST_F(PolygonTest, isPntInPolygonCheckPointsRestOnPolygonEdges)
{
for (std::size_t k(0); k<_pnts.size()-1; k++) {
for (double t(0.0); t<1.0; t+=0.001) {
@@ -94,20 +98,128 @@ TEST_F(IsPntInPolygonTest, CheckPointsRestOnPolygonEdges)
}
}
-TEST_F(IsPntInPolygonTest, CheckInnerPoints)
+TEST_F(PolygonTest, isPntInPolygonCheckInnerPoints)
{
ASSERT_TRUE(_polygon->isPntInPolygon(Point(1.0,1.0,0.0)));
ASSERT_TRUE(_polygon->isPntInPolygon(Point(0.5,1.0,0.0)));
}
-TEST_F(IsPntInPolygonTest, CheckOuterPoints)
+TEST_F(PolygonTest, isPntInPolygonCheckOuterPoints)
{
- ASSERT_FALSE(_polygon->isPntInPolygon(Point(1.0-std::numeric_limits<float>::epsilon(),0.0,0.0)));
- ASSERT_FALSE(_polygon->isPntInPolygon(Point(1.0+std::numeric_limits<float>::epsilon(),0.0,0.0)));
- ASSERT_FALSE(_polygon->isPntInPolygon(Point(0.0-std::numeric_limits<float>::epsilon(),1.0,0.0)));
- ASSERT_FALSE(_polygon->isPntInPolygon(Point(0.0-std::numeric_limits<float>::epsilon(),2.0,0.0)));
- ASSERT_FALSE(_polygon->isPntInPolygon(Point(1.0-std::numeric_limits<float>::epsilon(),3.0,0.0)));
- ASSERT_FALSE(_polygon->isPntInPolygon(Point(1.0+std::numeric_limits<float>::epsilon(),3.0,0.0)));
+ ASSERT_FALSE(_polygon->isPntInPolygon(Point(0.0-std::numeric_limits<float>::epsilon(),0.0,0.0)));
+ ASSERT_FALSE(_polygon->isPntInPolygon(Point(-2.0-std::numeric_limits<float>::epsilon(),2.0,0.0)));
+ ASSERT_FALSE(_polygon->isPntInPolygon(Point(-2.0-std::numeric_limits<float>::epsilon(),4.0,0.0)));
+ ASSERT_FALSE(_polygon->isPntInPolygon(Point(-1.0, 2.0+std::numeric_limits<float>::epsilon(),0.0)));
+ ASSERT_FALSE(_polygon->isPntInPolygon(Point(0.0-std::numeric_limits<float>::epsilon(),4.0,0.0)));
+ ASSERT_FALSE(_polygon->isPntInPolygon(Point(1.0,2.0+std::numeric_limits<float>::epsilon(),0.0)));
+ ASSERT_FALSE(_polygon->isPntInPolygon(Point(2.0-std::numeric_limits<float>::epsilon(),4.0,0.0)));
ASSERT_FALSE(_polygon->isPntInPolygon(Point(2.0+std::numeric_limits<float>::epsilon(),2.0,0.0)));
- ASSERT_FALSE(_polygon->isPntInPolygon(Point(2.0+std::numeric_limits<float>::epsilon(),1.0,0.0)));
+}
+
+/**
+ * 2 4 6
+ * |\ / \ /|
+ * | \ / \ / |
+ * 1 3 5 7
+ * \ /
+ * \ /
+ * \ /
+ * \ /
+ * \ /
+ * 0
+ * 0 = (0,0), 1=(-2,2), 2=(-2,4), 3=(-1,2), 4=(0,4), 5=(1,2), 6=(2,4), 7=(2,2)
+ */
+TEST_F(PolygonTest, containsSegment)
+{
+ // test segment (2,6)
+ GeoLib::Point &a(*(const_cast<GeoLib::Point*>(_polygon->getPoint(2))));
+ GeoLib::Point &b(*(const_cast<GeoLib::Point*>(_polygon->getPoint(6))));
+ ASSERT_FALSE(_polygon->containsSegment(a, b));
+
+ // test segments of polygon
+ for (std::size_t k(0); k<_polygon->getNumberOfPoints()-1;k++) {
+ a = *(const_cast<GeoLib::Point*>(_polygon->getPoint(k)));
+ b = *(const_cast<GeoLib::Point*>(_polygon->getPoint(k+1)));
+ EXPECT_TRUE(_polygon->containsSegment(a, b));
+ }
+
+ // 01
+ a = *(const_cast<GeoLib::Point*>(_polygon->getPoint(0)));
+ b = *(const_cast<GeoLib::Point*>(_polygon->getPoint(1)));
+ ASSERT_TRUE(_polygon->containsSegment(a, b));
+
+ // 12
+ a = *(const_cast<GeoLib::Point*>(_polygon->getPoint(1)));
+ b = *(const_cast<GeoLib::Point*>(_polygon->getPoint(2)));
+ ASSERT_TRUE(_polygon->containsSegment(a, b));
+
+ // 23
+ a = *(const_cast<GeoLib::Point*>(_polygon->getPoint(2)));
+ b = *(const_cast<GeoLib::Point*>(_polygon->getPoint(3)));
+ ASSERT_TRUE(_polygon->containsSegment(a, b));
+
+ // 34
+ a = *(const_cast<GeoLib::Point*>(_polygon->getPoint(3)));
+ b = *(const_cast<GeoLib::Point*>(_polygon->getPoint(4)));
+ ASSERT_TRUE(_polygon->containsSegment(a, b));
+
+ // 45
+ a = *(const_cast<GeoLib::Point*>(_polygon->getPoint(4)));
+ b = *(const_cast<GeoLib::Point*>(_polygon->getPoint(5)));
+ ASSERT_TRUE(_polygon->containsSegment(a, b));
+
+ // 56
+ a = *(const_cast<GeoLib::Point*>(_polygon->getPoint(5)));
+ b = *(const_cast<GeoLib::Point*>(_polygon->getPoint(6)));
+ ASSERT_TRUE(_polygon->containsSegment(a, b));
+
+ // 67
+ a = *(const_cast<GeoLib::Point*>(_polygon->getPoint(6)));
+ b = *(const_cast<GeoLib::Point*>(_polygon->getPoint(7)));
+ ASSERT_TRUE(_polygon->containsSegment(a, b));
+
+ // 70
+ a = *(const_cast<GeoLib::Point*>(_polygon->getPoint(7)));
+ b = *(const_cast<GeoLib::Point*>(_polygon->getPoint(0)));
+ ASSERT_TRUE(_polygon->containsSegment(a, b));
+
+ // test segment (3,5)
+ a = *(const_cast<GeoLib::Point*>(_polygon->getPoint(3)));
+ b = *(const_cast<GeoLib::Point*>(_polygon->getPoint(5)));
+ ASSERT_TRUE(_polygon->containsSegment(a, b));
+
+ // test segment (1,7)
+ a = *(const_cast<GeoLib::Point*>(_polygon->getPoint(1)));
+ b = *(const_cast<GeoLib::Point*>(_polygon->getPoint(7)));
+ ASSERT_TRUE(_polygon->containsSegment(a, b));
+
+ // test segment (1,4)
+ a = *(const_cast<GeoLib::Point*>(_polygon->getPoint(1)));
+ b = *(const_cast<GeoLib::Point*>(_polygon->getPoint(4)));
+ ASSERT_FALSE(_polygon->containsSegment(a, b));
+
+}
+
+TEST_F(PolygonTest, isPolylineInPolygon)
+{
+ // create a test polyline
+ std::vector<GeoLib::Point*> pnts;
+ pnts.push_back(new GeoLib::Point(-2.0,4.0,0.0)); // 2
+ pnts.push_back(new GeoLib::Point( 2.0,4.0,0.0)); // 6
+ GeoLib::Polyline outer_ply(pnts);
+ outer_ply.addPoint(0);
+ outer_ply.addPoint(1);
+ ASSERT_FALSE(_polygon->isPolylineInPolygon(outer_ply));
+ for (std::size_t k(0); k<pnts.size(); k++)
+ delete pnts[k];
+ pnts.clear();
+
+ pnts.push_back(new GeoLib::Point(-1.0,2.0,0.0)); // 3
+ pnts.push_back(new GeoLib::Point( 1.0,2.0,0.0)); // 5
+ GeoLib::Polyline inner_ply(pnts);
+ inner_ply.addPoint(0);
+ inner_ply.addPoint(1);
+ ASSERT_TRUE(_polygon->isPolylineInPolygon(inner_ply));
+ for (std::size_t k(0); k<pnts.size(); k++)
+ delete pnts[k];
}
diff --git a/Tests/GeoLib/TestSimplePolygonTree.cpp b/Tests/GeoLib/TestSimplePolygonTree.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..aa31f022cfff9731da849eb4bd73d39c00ccc0dc
--- /dev/null
+++ b/Tests/GeoLib/TestSimplePolygonTree.cpp
@@ -0,0 +1,177 @@
+/**
+ * @file TestSimplePolygonTree.cpp
+ * @author
+ * @date Apr 02, 2014
+ * @brief Test functionality of class SimplePolygonTree.
+ *
+ * @copyright
+ * Copyright (c) 2014, OpenGeoSys Community (http://www.opengeosys.org)
+ * Distributed under a Modified BSD License.
+ * See accompanying file LICENSE.txt or
+ * http://www.opengeosys.org/LICENSE.txt
+ */
+
+#include <memory>
+
+#include "gtest/gtest.h"
+
+#include "Point.h"
+#include "Polygon.h"
+#include "SimplePolygonTree.h"
+
+using namespace GeoLib;
+
+/**
+ * 2
+ * / \
+ * / \
+ * / \
+ * / \
+ * 3---------1
+ * \ 6--5 /
+ * \ \ / /
+ * \ 4 /
+ * \ /
+ * 0
+ *
+ * Polygons: P0: 3,1,2,3 / P1: 0,1,3,0 / P2: 0,1,2,3,0
+ * P3: 4,5,6,4
+ */
+
+
+class CreatePolygonTreesTest : public testing::Test
+{
+public:
+ CreatePolygonTreesTest() :
+ _p0(nullptr), _p1(nullptr), _p2(nullptr), _p3(nullptr)
+ {
+ // create points and construct polygon
+ _pnts.push_back(new Point(0.0,-1.0,0.0));
+ _pnts.push_back(new Point(1.0,0.0,0.0));
+ _pnts.push_back(new Point(0.0,1.0,0.0));
+ _pnts.push_back(new Point(-1.0,0.0,0.0));
+
+ _pnts.push_back(new Point(-0.9,0.0,0.0));
+ _pnts.push_back(new Point(0.75,-0.1,0.0));
+ _pnts.push_back(new Point(-0.75,-0.1,0.0));
+
+ // create closed polylines
+ Polyline ply0(_pnts);
+ ply0.addPoint(3);
+ ply0.addPoint(1);
+ ply0.addPoint(2);
+ ply0.addPoint(3);
+
+ Polyline ply1(_pnts);
+ ply1.addPoint(0);
+ ply1.addPoint(1);
+ ply1.addPoint(3);
+ ply1.addPoint(0);
+
+ Polyline ply2(_pnts);
+ ply2.addPoint(0);
+ ply2.addPoint(1);
+ ply2.addPoint(2);
+ ply2.addPoint(3);
+ ply2.addPoint(0);
+
+ Polyline ply3(_pnts);
+ ply3.addPoint(4);
+ ply3.addPoint(5);
+ ply3.addPoint(6);
+ ply3.addPoint(4);
+
+ // create polygons
+ _p0 = new Polygon(ply0);
+ _p1 = new Polygon(ply1);
+ _p2 = new Polygon(ply2);
+ _p3 = new Polygon(ply3);
+ }
+
+ ~CreatePolygonTreesTest()
+ {
+ delete _p0;
+ delete _p1;
+ delete _p2;
+ delete _p3;
+ for (std::size_t k(0); k<_pnts.size(); k++)
+ delete _pnts[k];
+ }
+
+protected:
+ std::vector<Point*> _pnts;
+ Polygon *_p0;
+ Polygon *_p1;
+ Polygon *_p2;
+ Polygon *_p3;
+};
+
+TEST_F(CreatePolygonTreesTest, P0AndP1)
+{
+ SimplePolygonTree *pt0(new SimplePolygonTree(_p0, nullptr));
+ SimplePolygonTree *pt1(new SimplePolygonTree(_p1, nullptr));
+
+ std::list<SimplePolygonTree*> pt_list;
+ pt_list.push_back(pt0);
+ pt_list.push_back(pt1);
+ ASSERT_EQ(2u, pt_list.size());
+
+ createPolygonTrees(pt_list);
+ ASSERT_EQ(2u, pt_list.size());
+ std::for_each(pt_list.begin(), pt_list.end(), std::default_delete<GeoLib::SimplePolygonTree>());
+}
+
+TEST_F(CreatePolygonTreesTest, P0AndP1AndP2)
+{
+ SimplePolygonTree *pt0(new SimplePolygonTree(_p0, nullptr));
+ SimplePolygonTree *pt1(new SimplePolygonTree(_p1, nullptr));
+ SimplePolygonTree *pt2(new SimplePolygonTree(_p2, nullptr));
+
+ std::list<SimplePolygonTree*> pt_list;
+ pt_list.push_back(pt0);
+ pt_list.push_back(pt1);
+ pt_list.push_back(pt2);
+ ASSERT_EQ(3u, pt_list.size());
+
+ ASSERT_FALSE(_p0->isPolylineInPolygon(*_p1));
+ ASSERT_FALSE(_p0->isPolylineInPolygon(*_p2));
+ ASSERT_FALSE(_p1->isPolylineInPolygon(*_p0));
+ ASSERT_FALSE(_p1->isPolylineInPolygon(*_p2));
+ ASSERT_TRUE(_p2->isPolylineInPolygon(*_p0));
+ ASSERT_TRUE(_p2->isPolylineInPolygon(*_p1));
+
+ createPolygonTrees(pt_list);
+ ASSERT_EQ(1u, pt_list.size());
+ ASSERT_EQ(2u, (*(pt_list.begin()))->getNChilds());
+ std::for_each(pt_list.begin(), pt_list.end(), std::default_delete<GeoLib::SimplePolygonTree>());
+}
+
+TEST_F(CreatePolygonTreesTest, P0AndP1AndP2AndP3)
+{
+ SimplePolygonTree *pt0(new SimplePolygonTree(_p0, nullptr));
+ SimplePolygonTree *pt1(new SimplePolygonTree(_p1, nullptr));
+ SimplePolygonTree *pt2(new SimplePolygonTree(_p2, nullptr));
+ SimplePolygonTree *pt3(new SimplePolygonTree(_p3, nullptr));
+
+ std::list<SimplePolygonTree*> pt_list;
+ pt_list.push_back(pt0);
+ pt_list.push_back(pt1);
+ pt_list.push_back(pt2);
+ pt_list.push_back(pt3);
+ ASSERT_EQ(4u, pt_list.size());
+
+ ASSERT_FALSE(_p0->isPolylineInPolygon(*_p1));
+ ASSERT_FALSE(_p0->isPolylineInPolygon(*_p2));
+ ASSERT_FALSE(_p1->isPolylineInPolygon(*_p0));
+ ASSERT_FALSE(_p1->isPolylineInPolygon(*_p2));
+ ASSERT_TRUE(_p2->isPolylineInPolygon(*_p0));
+ ASSERT_TRUE(_p2->isPolylineInPolygon(*_p1));
+ ASSERT_TRUE(_p2->isPolylineInPolygon(*_p3));
+ ASSERT_TRUE(_p1->isPolylineInPolygon(*_p3));
+
+ createPolygonTrees(pt_list);
+ ASSERT_EQ(1u, pt_list.size());
+ ASSERT_EQ(2u, (*(pt_list.begin()))->getNChilds());
+ std::for_each(pt_list.begin(), pt_list.end(), std::default_delete<GeoLib::SimplePolygonTree>());
+}
+
diff --git a/Tests/MathLib/TestGlobalMatrixInterface.cpp b/Tests/MathLib/TestGlobalMatrixInterface.cpp
index 04e1fe8637506b36b252e0e83d1082d67a84f47b..d11f17493692d6f85bca41a19799d52f9f4a33d8 100644
--- a/Tests/MathLib/TestGlobalMatrixInterface.cpp
+++ b/Tests/MathLib/TestGlobalMatrixInterface.cpp
@@ -178,7 +178,7 @@ TEST(Math, CheckInterface_LisMatrix)
#endif
#ifdef USE_PETSC // or MPI
-TEST(Math, CheckInterface_PETScMatrix_Local_Size)
+TEST(MPITest_Math, CheckInterface_PETScMatrix_Local_Size)
{
MathLib::PETScMatrixOption opt;
opt.d_nz = 2;
@@ -193,7 +193,7 @@ TEST(Math, CheckInterface_PETScMatrix_Local_Size)
checkGlobalMatrixInterfaceMPI(A, x);
}
-TEST(Math, CheckInterface_PETScMatrix_Global_Size)
+TEST(MPITest_Math, CheckInterface_PETScMatrix_Global_Size)
{
MathLib::PETScMatrixOption opt;
opt.d_nz = 2;
@@ -206,7 +206,7 @@ TEST(Math, CheckInterface_PETScMatrix_Global_Size)
}
// Test rectangular matrix
-TEST(Math, CheckInterface_PETSc_Rectangular_Matrix_Local_Size)
+TEST(MPITest_Math, CheckInterface_PETSc_Rectangular_Matrix_Local_Size)
{
MathLib::PETScMatrixOption opt;
opt.d_nz = 3;
@@ -221,7 +221,7 @@ TEST(Math, CheckInterface_PETSc_Rectangular_Matrix_Local_Size)
checkGlobalRectangularMatrixInterfaceMPI(A, x);
}
-TEST(Math, CheckInterface_PETSc_Rectangular_Matrix_Global_Size)
+TEST(MPITest_Math, CheckInterface_PETSc_Rectangular_Matrix_Global_Size)
{
MathLib::PETScMatrixOption opt;
opt.d_nz = 3;
diff --git a/Tests/MathLib/TestGlobalVectorInterface.cpp b/Tests/MathLib/TestGlobalVectorInterface.cpp
index bfcf01796a44a2527f192a07874b553f8f5451e4..01e21e28c164af2176944cef6cde7dc75052d39f 100644
--- a/Tests/MathLib/TestGlobalVectorInterface.cpp
+++ b/Tests/MathLib/TestGlobalVectorInterface.cpp
@@ -204,7 +204,7 @@ TEST(Math, CheckInterface_LisVector)
//--------------------------------------------
#ifdef USE_PETSC
-TEST(Math, CheckInterface_PETScVector)
+TEST(MPITest_Math, CheckInterface_PETScVector)
{
checkGlobalVectorInterfaceMPI<MathLib::PETScVector >();
}
diff --git a/Tests/MathLib/TestLinearSolver.cpp b/Tests/MathLib/TestLinearSolver.cpp
index 6c46e62fe383adcbf65a7ee3b5886d4a0d0ba4e3..a8786535c962ef177cf2e09aa261901de1344c0e 100644
--- a/Tests/MathLib/TestLinearSolver.cpp
+++ b/Tests/MathLib/TestLinearSolver.cpp
@@ -1,7 +1,8 @@
/**
* \file
* \author Norihiro Watanabe
- * \date 2013-04-16
+ * \author Wenqing Wang
+ * \date 2013-04-16, 2014-04
* \brief Implementation tests.
*
* \copyright
@@ -21,11 +22,19 @@
#include "MathLib/LinAlg/Dense/DenseTools.h"
#include "MathLib/LinAlg/FinalizeMatrixAssembly.h"
#include "MathLib/LinAlg/Solvers/GaussAlgorithm.h"
+
#ifdef USE_LIS
#include "MathLib/LinAlg/Lis/LisLinearSolver.h"
#include "MathLib/LinAlg/Lis/LisTools.h"
#endif
+#ifdef USE_PETSC
+#include "MathLib/LinAlg/PETSc/PETScMatrix.h"
+#include "MathLib/LinAlg/PETSc/PETScVector.h"
+#include "MathLib/LinAlg/PETSc/PETScLinearSolver.h"
+#include "MathLib/LinAlg/PETSc/PETScTools.h"
+#endif
+
#include "../TestTools.h"
namespace
@@ -34,7 +43,8 @@ namespace
template<class T_Mat>
void setMatrix9x9(T_Mat &mat)
{
- double d_mat[] = {
+ double d_mat[] =
+ {
6.66667e-012, -1.66667e-012, 0, -1.66667e-012, -3.33333e-012, 0, 0, 0, 0,
-1.66667e-012, 1.33333e-011, -1.66667e-012, -3.33333e-012, -3.33333e-012, -3.33333e-012, 0, 0, 0,
0, -1.66667e-012, 6.66667e-012, 0, -3.33333e-012, -1.66667e-012, 0, 0, 0,
@@ -45,9 +55,9 @@ void setMatrix9x9(T_Mat &mat)
0, 0, 0, -3.33333e-012, -3.33333e-012, -3.33333e-012, -1.66667e-012, 1.33333e-011, -1.66667e-012,
0, 0, 0, 0, -3.33333e-012, -1.66667e-012, 0, -1.66667e-012, 6.66667e-012
};
- for (unsigned i = 0; i < 9; i++)
- for (unsigned j = 0; j < 9; j++)
- mat.setValue(i, j, d_mat[i*9+j]);
+ for (unsigned i = 0; i < 9; i++)
+ for (unsigned j = 0; j < 9; j++)
+ mat.setValue(i, j, d_mat[i*9+j]);
}
struct Example1
@@ -59,7 +69,7 @@ struct Example1
double* exH;
Example1()
- : mat(dim_eqs, dim_eqs), exH(new double[dim_eqs])
+ : mat(dim_eqs, dim_eqs), exH(new double[dim_eqs])
{
setMatrix9x9(mat);
std::size_t int_dirichlet_bc_id[] = {2,5,8,0,3,6};
@@ -67,7 +77,8 @@ struct Example1
vec_dirichlet_bc_value.resize(6);
std::fill(vec_dirichlet_bc_value.begin(), vec_dirichlet_bc_value.begin()+3, .0);
std::fill(vec_dirichlet_bc_value.begin()+3, vec_dirichlet_bc_value.end(), 1.0);
- for (std::size_t i=0; i<9; i++) {
+ for (std::size_t i=0; i<9; i++)
+ {
if (i%3==0) exH[i] = 1.0;
if (i%3==1) exH[i] = 0.5;
if (i%3==2) exH[i] = 0.;
@@ -87,8 +98,10 @@ void checkLinearSolverInterface(T_MATRIX &A, boost::property_tree::ptree &ls_opt
// set a coefficient matrix
A.setZero();
- for (size_t i=0; i<ex1.dim_eqs; i++) {
- for (size_t j=0; j<ex1.dim_eqs; j++) {
+ for (size_t i=0; i<ex1.dim_eqs; i++)
+ {
+ for (size_t j=0; j<ex1.dim_eqs; j++)
+ {
double v = ex1.mat(i, j);
if (v!=.0)
A.add(i, j, v);
@@ -112,6 +125,71 @@ void checkLinearSolverInterface(T_MATRIX &A, boost::property_tree::ptree &ls_opt
}
+#ifdef USE_PETSC
+template <class T_MATRIX, class T_VECTOR, class T_LINEAR_SOVLER>
+void checkLinearSolverInterface(T_MATRIX &A, T_VECTOR &b, const std::string &prefix_name)
+{
+ int mrank;
+ MPI_Comm_rank(PETSC_COMM_WORLD, &mrank);
+ // Add entries
+ MathLib::DenseMatrix<double> loc_m(2, 2);
+ loc_m(0, 0) = 1. + mrank;
+ loc_m(0, 1) = 2. + mrank;
+ loc_m(1, 0) = 3. + mrank;
+ loc_m(1, 1) = 4. + mrank;
+
+ std::vector<int> row_pos(2);
+ std::vector<int> col_pos(2);
+ row_pos[0] = 2 * mrank;
+ row_pos[1] = 2 * mrank + 1;
+ col_pos[0] = row_pos[0];
+ col_pos[1] = row_pos[1];
+
+ A.add(row_pos, col_pos, loc_m);
+
+ MathLib::finalizeMatrixAssembly(A);
+
+ const bool deep_copy = false;
+ T_VECTOR x(b, deep_copy);
+
+ std::vector<double> local_vec(2);
+ local_vec[0] = mrank+1;
+ local_vec[1] = 2. * (mrank+1);
+ x.set(row_pos, local_vec);
+
+ double x0[6];
+ double x1[6];
+ x.getGlobalVector(x0);
+
+ A.multiply(x, b);
+
+ // apply BC
+ std::vector<int> bc_id; // Type must be int to match Petsc_Int
+ std::vector<double> bc_value;
+
+ if(mrank == 1)
+ {
+ bc_id.resize(1);
+ bc_value.resize(1);
+ bc_id[0] = 2 * mrank;
+ bc_value[0] = mrank+1;
+ }
+
+ MathLib::applyKnownSolution(A, b, x, bc_id, bc_value);
+
+ MathLib::finalizeMatrixAssembly(A);
+
+ // solve
+ T_LINEAR_SOVLER ls(A, prefix_name);
+ EXPECT_TRUE(ls.solve(b, x));
+
+ EXPECT_GT(ls.getNumberOfIterations(), 0u);
+
+ x.getGlobalVector(x1);
+ ASSERT_ARRAY_NEAR(x0, x1, 6, 1e-5);
+}
+#endif
+
} // end namespace
TEST(MathLib, CheckInterface_GaussAlgorithm)
@@ -142,3 +220,80 @@ TEST(Math, CheckInterface_Lis)
}
#endif
+#ifdef USE_PETSC
+TEST(MPITest_Math, CheckInterface_PETSc_Linear_Solver_basic)
+{
+ MathLib::PETScMatrixOption opt;
+ opt.d_nz = 2;
+ opt.o_nz = 0;
+ opt.is_global_size = false;
+ opt.n_local_cols = 2;
+ MathLib::PETScMatrix A(2, opt);
+
+ const bool is_gloabal_size = false;
+ MathLib::PETScVector b(2, is_gloabal_size);
+
+ PetscOptionsSetValue("-ptest1_ksp_type", "bcgs");
+
+ PetscOptionsSetValue("-ptest1_ksp_rtol", "1.e-8");
+ PetscOptionsSetValue("-ptest1_ksp_atol", "1.e-50");
+ PetscOptionsSetValue("-ptest1_ksp_max_it", "1000");
+
+ PetscOptionsSetValue("-ptest1_pc_type", "bjacobi");
+
+ checkLinearSolverInterface<MathLib::PETScMatrix, MathLib::PETScVector,
+ MathLib::PETScLinearSolver>(A, b, "ptest1_");
+}
+
+TEST(MPITest_Math, CheckInterface_PETSc_Linear_Solver_chebyshev_sor)
+{
+ MathLib::PETScMatrixOption opt;
+ opt.d_nz = 2;
+ opt.o_nz = 0;
+ opt.is_global_size = false;
+ opt.n_local_cols = 2;
+ MathLib::PETScMatrix A(2, opt);
+
+ const bool is_gloabal_size = false;
+ MathLib::PETScVector b(2, is_gloabal_size);
+
+ PetscOptionsSetValue("-ptest2_ksp_type", "chebyshev");
+
+ PetscOptionsSetValue("-ptest2_ksp_rtol", "1.e-8");
+ PetscOptionsSetValue("-ptest2_ksp_atol", "1.e-50");
+ PetscOptionsSetValue("-ptest2_ksp_max_it", "1000");
+
+ PetscOptionsSetValue("-ptest2_pc_type", "sor");
+
+ checkLinearSolverInterface<MathLib::PETScMatrix, MathLib::PETScVector,
+ MathLib::PETScLinearSolver>(A, b, "ptest2_");
+}
+
+TEST(MPITest_Math, CheckInterface_PETSc_Linear_Solver_gmres_amg)
+{
+ MathLib::PETScMatrixOption opt;
+ opt.d_nz = 2;
+ opt.o_nz = 0;
+ opt.is_global_size = false;
+ opt.n_local_cols = 2;
+ MathLib::PETScMatrix A(2, opt);
+
+ const bool is_gloabal_size = false;
+ MathLib::PETScVector b(2, is_gloabal_size);
+
+ PetscOptionsSetValue("-ptest3_ksp_type", "gmres");
+ PetscOptionsSetValue("-ptest3_ksp_rtol", "1.e-8");
+ PetscOptionsSetValue("-ptest3_ksp_gmres_restart", "20");
+ PetscOptionsSetValue("-ptest3_ksp_gmres_classicalgramschmidt", "");
+
+ PetscOptionsSetValue("-ptest3_pc_type", "gamg");
+ PetscOptionsSetValue("-ptest3_pc_gamg_type", "agg");
+ PetscOptionsSetValue("-ptest3_pc_gamg_agg_nsmooths", "2");
+
+ checkLinearSolverInterface<MathLib::PETScMatrix, MathLib::PETScVector,
+ MathLib::PETScLinearSolver>(A, b, "ptest3_");
+}
+
+#endif
+
+
diff --git a/Utils/MeshEdit/CMakeLists.txt b/Utils/MeshEdit/CMakeLists.txt
index 908626e2986f94bb6adec1db1153f353395089cd..12101c8b13a627e6a12300ffe4b30a63a90a81be 100644
--- a/Utils/MeshEdit/CMakeLists.txt
+++ b/Utils/MeshEdit/CMakeLists.txt
@@ -57,3 +57,14 @@ TARGET_LINK_LIBRARIES( MoveMesh
)
SET_TARGET_PROPERTIES(MoveMesh PROPERTIES FOLDER Utilities)
+ADD_EXECUTABLE( appendLinesAlongPolyline appendLinesAlongPolyline.cpp )
+TARGET_LINK_LIBRARIES( appendLinesAlongPolyline
+ BaseLib
+ FileIO
+ MathLib
+ MeshLib
+ MeshGeoToolsLib
+)
+
+SET_TARGET_PROPERTIES(appendLinesAlongPolyline
+ PROPERTIES FOLDER Utilities)
diff --git a/Utils/MeshEdit/appendLinesAlongPolyline.cpp b/Utils/MeshEdit/appendLinesAlongPolyline.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..d4a86db5dd3ac5f04c31ce0e38c2a53e99dc9798
--- /dev/null
+++ b/Utils/MeshEdit/appendLinesAlongPolyline.cpp
@@ -0,0 +1,103 @@
+/**
+ * @copyright
+ * Copyright (c) 2013, OpenGeoSys Community (http://www.opengeosys.org)
+ * Distributed under a Modified BSD License.
+ * See accompanying file LICENSE.txt or
+ * http://www.opengeosys.org/LICENSE.txt
+ */
+
+// TCLAP
+#include "tclap/CmdLine.h"
+
+// ThirdParty/logog
+#include "logog/include/logog.hpp"
+
+// BaseLib
+#include "LogogSimpleFormatter.h"
+#include "FileTools.h"
+
+// GeoLib
+#include "GEOObjects.h"
+#include "PolylineVec.h"
+
+// FileIO
+#include "Legacy/MeshIO.h"
+#include "readMeshFromFile.h"
+#include "XmlIO/Boost/BoostXmlGmlInterface.h"
+
+// MeshLib
+#include "Mesh.h"
+
+// MeshGeoToolsLib
+#include "MeshGeoToolsLib/AppendLinesAlongPolyline.h"
+
+
+int main (int argc, char* argv[])
+{
+ LOGOG_INITIALIZE();
+ logog::Cout* logog_cout (new logog::Cout);
+ BaseLib::LogogSimpleFormatter *custom_format (new BaseLib::LogogSimpleFormatter);
+ logog_cout->SetFormatter(*custom_format);
+
+ TCLAP::CmdLine cmd("Append line elements into a mesh.", ' ', "0.1");
+ TCLAP::ValueArg<std::string> mesh_in("i", "mesh-input-file",
+ "the name of the file containing the input mesh", true,
+ "", "file name of input mesh");
+ cmd.add(mesh_in);
+ TCLAP::ValueArg<std::string> mesh_out("o", "mesh-output-file",
+ "the name of the file the mesh will be written to", true,
+ "", "file name of output mesh");
+ cmd.add(mesh_out);
+ TCLAP::ValueArg<std::string> geoFileArg("g", "geo-file",
+ "the name of the geometry file which contains polylines", true, "", "the name of the geometry file");
+ cmd.add(geoFileArg);
+
+ // parse arguments
+ cmd.parse(argc, argv);
+
+ // read GEO objects
+ GeoLib::GEOObjects geo_objs;
+ FileIO::BoostXmlGmlInterface xml(geo_objs);
+ xml.readFile(geoFileArg.getValue());
+
+ std::vector<std::string> geo_names;
+ geo_objs.getGeometryNames (geo_names);
+ if (geo_names.empty ())
+ {
+ std::cout << "no geometries found" << std::endl;
+ return -1;
+ }
+ const GeoLib::PolylineVec* ply_vec (geo_objs.getPolylineVecObj(geo_names[0]));
+ if (!ply_vec)
+ {
+ std::cout << "could not found polylines" << std::endl;
+ return -1;
+ }
+
+ // read a mesh
+ MeshLib::Mesh const*const mesh (FileIO::readMeshFromFile(mesh_in.getValue()));
+ if (!mesh)
+ {
+ ERR("Mesh file %s not found", mesh_in.getValue().c_str());
+ return 1;
+ }
+ INFO("Mesh read: %d nodes, %d elements.", mesh->getNNodes(), mesh->getNElements());
+
+ // add line elements
+ MeshLib::Mesh* new_mesh = MeshGeoToolsLib::appendLinesAlongPolylines(*mesh, *ply_vec);
+ INFO("Mesh created: %d nodes, %d elements.", new_mesh->getNNodes(), new_mesh->getNElements());
+
+ // write into a file
+ FileIO::Legacy::MeshIO meshIO;
+ meshIO.setMesh(new_mesh);
+ meshIO.writeToFile(mesh_out.getValue());
+
+ delete new_mesh;
+
+ delete custom_format;
+ delete logog_cout;
+ LOGOG_SHUTDOWN();
+
+ return 1;
+}
+
diff --git a/Utils/SimpleMeshCreation/CMakeLists.txt b/Utils/SimpleMeshCreation/CMakeLists.txt
index 3127d017eca514b5d221f67b41af38774d79d2ee..b2ec276f40c8c2ec5e872c9237d413843509a11d 100644
--- a/Utils/SimpleMeshCreation/CMakeLists.txt
+++ b/Utils/SimpleMeshCreation/CMakeLists.txt
@@ -1,15 +1,16 @@
-IF(TARGET VtkVis)
- INCLUDE_DIRECTORIES(
- ${CMAKE_SOURCE_DIR}
- ${CMAKE_SOURCE_DIR}/BaseLib
- ${CMAKE_SOURCE_DIR}/FileIO
- ${CMAKE_SOURCE_DIR}/FileIO/Legacy
- ${CMAKE_SOURCE_DIR}/GeoLib
- ${CMAKE_SOURCE_DIR}/MathLib
- ${CMAKE_SOURCE_DIR}/MeshLib
- ${CMAKE_SOURCE_DIR}/Gui/VtkVis
- )
+INCLUDE_DIRECTORIES(
+ ${CMAKE_SOURCE_DIR}
+ ${CMAKE_SOURCE_DIR}/BaseLib
+ ${CMAKE_SOURCE_DIR}/FileIO
+ ${CMAKE_SOURCE_DIR}/FileIO/Legacy
+ ${CMAKE_SOURCE_DIR}/GeoLib
+ ${CMAKE_SOURCE_DIR}/MathLib
+ ${CMAKE_SOURCE_DIR}/MeshLib
+ ${CMAKE_SOURCE_DIR}/Gui/VtkVis
+)
+
+IF(TARGET VtkVis)
ADD_EXECUTABLE( generateStructuredQuadMesh generateStructuredQuadMesh.cpp )
SET_TARGET_PROPERTIES( generateStructuredQuadMesh PROPERTIES FOLDER Utils)
@@ -35,3 +36,14 @@ IF(TARGET VtkVis)
)
ENDIF() # VtkVis-target is existing
+
+
+ADD_EXECUTABLE( generateStructuredMesh generateStructuredMesh.cpp )
+TARGET_LINK_LIBRARIES( generateStructuredMesh
+ BaseLib
+ FileIO
+ MathLib
+ MeshLib
+)
+SET_TARGET_PROPERTIES(generateStructuredMesh PROPERTIES FOLDER Utilities)
+
diff --git a/Utils/SimpleMeshCreation/generateStructuredMesh.cpp b/Utils/SimpleMeshCreation/generateStructuredMesh.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..f25e05b37c508991485d203c27384b5e6227f03f
--- /dev/null
+++ b/Utils/SimpleMeshCreation/generateStructuredMesh.cpp
@@ -0,0 +1,97 @@
+/**
+ * @copyright
+ * Copyright (c) 2013, OpenGeoSys Community (http://www.opengeosys.org)
+ * Distributed under a Modified BSD License.
+ * See accompanying file LICENSE.txt or
+ * http://www.opengeosys.org/LICENSE.txt
+ */
+
+// TCLAP
+#include "tclap/CmdLine.h"
+
+// ThirdParty/logog
+#include "logog/include/logog.hpp"
+
+// BaseLib
+#include "LogogSimpleFormatter.h"
+
+// FileIO
+#include "Legacy/MeshIO.h"
+
+// MeshLib
+#include "Mesh.h"
+#include "Node.h"
+#include "Elements/Element.h"
+#include "MeshEnums.h"
+#include "MeshGenerators/MeshGenerator.h"
+
+
+int main (int argc, char* argv[])
+{
+ LOGOG_INITIALIZE();
+ logog::Cout* logog_cout (new logog::Cout);
+ BaseLib::LogogSimpleFormatter *custom_format (new BaseLib::LogogSimpleFormatter);
+ logog_cout->SetFormatter(*custom_format);
+
+ TCLAP::CmdLine cmd("Generate a structured mesh.", ' ', "0.1");
+ TCLAP::ValueArg<std::string> mesh_out("o", "mesh-output-file",
+ "the name of the file the mesh will be written to", true,
+ "", "file name of output mesh");
+ cmd.add(mesh_out);
+ TCLAP::ValueArg<std::string> eleTypeArg("e", "element-type",
+ "element type to be created", true, "line", "element type");
+ cmd.add(eleTypeArg);
+ TCLAP::ValueArg<double> lengthArg("l", "length",
+ "length of a domain", true, 10.0, "length of a domain");
+ cmd.add(lengthArg);
+ TCLAP::ValueArg<unsigned> nsubdivArg("n", "nr-subdivision",
+ "the number of subdivision", true, 10, "the number of subdivision");
+ cmd.add(nsubdivArg);
+
+ // parse arguments
+ cmd.parse(argc, argv);
+ const std::string eleTypeName(eleTypeArg.getValue());
+ const MeshElemType eleType = String2MeshElemType(eleTypeName);
+ const double length = lengthArg.getValue();
+ const unsigned n_subdivision = nsubdivArg.getValue();
+
+ // generate a mesh
+ MeshLib::Mesh* mesh = nullptr;
+ switch (eleType)
+ {
+ case MeshElemType::LINE:
+ mesh = MeshLib::MeshGenerator::generateLineMesh(length, n_subdivision);
+ break;
+ case MeshElemType::TRIANGLE:
+ mesh = MeshLib::MeshGenerator::generateRegularTriMesh(length, n_subdivision);
+ break;
+ case MeshElemType::QUAD:
+ mesh = MeshLib::MeshGenerator::generateRegularQuadMesh(length, n_subdivision);
+ break;
+ case MeshElemType::HEXAHEDRON:
+ mesh = MeshLib::MeshGenerator::generateRegularHexMesh(length, n_subdivision);
+ break;
+ default:
+ ERR("Given element type is not supported.");
+ break;
+ }
+
+ if (mesh)
+ {
+ INFO("Mesh created: %d nodes, %d elements.", mesh->getNNodes(), mesh->getNElements());
+
+ // write into a file
+ FileIO::Legacy::MeshIO meshIO;
+ meshIO.setMesh(mesh);
+ meshIO.writeToFile(mesh_out.getValue());
+
+ delete mesh;
+ }
+
+ delete custom_format;
+ delete logog_cout;
+ LOGOG_SHUTDOWN();
+
+ return 1;
+}
+