diff --git a/GeoLib/AnalyticalGeometry-impl.h b/GeoLib/AnalyticalGeometry-impl.h
index 76090a3675deb1ee9ad0b98cc659fe18215e14c0..cf5596402d65e07474a2d2628f06a7b696285612 100644
--- a/GeoLib/AnalyticalGeometry-impl.h
+++ b/GeoLib/AnalyticalGeometry-impl.h
@@ -70,9 +70,8 @@ std::pair<Eigen::Vector3d, double> getNewellPlane(
return std::make_pair(plane_normal, d);
}
-template<class T_MATRIX>
-void compute2DRotationMatrixToX(MathLib::Vector3 const& v,
- T_MATRIX & rot_mat)
+template <class T_MATRIX>
+void compute2DRotationMatrixToX(Eigen::Vector3d const& v, T_MATRIX& rot_mat)
{
const double cos_theta = v[0];
const double sin_theta = v[1];
@@ -84,35 +83,43 @@ void compute2DRotationMatrixToX(MathLib::Vector3 const& v,
}
template <class T_MATRIX>
-void compute3DRotationMatrixToX(MathLib::Vector3 const& v,
- T_MATRIX & rot_mat)
+void compute3DRotationMatrixToX(Eigen::Vector3d const& v, T_MATRIX& rot_mat)
{
+ double const eps = std::numeric_limits<double>::epsilon();
// a vector on the plane
- MathLib::Vector3 yy(0, 0, 0);
- if (std::abs(v[0]) > 0.0 && std::abs(v[1]) + std::abs(v[2]) < std::numeric_limits<double>::epsilon()) {
+ Eigen::Vector3d yy({0, 0, 0});
+ if (std::abs(v[0]) > 0.0 && std::abs(v[1]) + std::abs(v[2]) < eps)
+ {
yy[2] = 1.0;
- } else if (std::abs(v[1]) > 0.0 && std::abs(v[0]) + std::abs(v[2]) < std::numeric_limits<double>::epsilon()) {
+ }
+ else if (std::abs(v[1]) > 0.0 && std::abs(v[0]) + std::abs(v[2]) < eps)
+ {
yy[0] = 1.0;
- } else if (std::abs(v[2]) > 0.0 && std::abs(v[0]) + std::abs(v[1]) < std::numeric_limits<double>::epsilon()) {
+ }
+ else if (std::abs(v[2]) > 0.0 && std::abs(v[0]) + std::abs(v[1]) < eps)
+ {
yy[1] = 1.0;
- } else {
- for (unsigned i = 0; i < 3; i++) {
- if (std::abs(v[i]) > 0.0) {
+ }
+ else
+ {
+ for (unsigned i = 0; i < 3; i++)
+ {
+ if (std::abs(v[i]) > 0.0)
+ {
yy[i] = -v[i];
break;
}
}
}
// z"_vec
- MathLib::Vector3 zz(MathLib::crossProduct(v, yy));
- zz.normalize();
+ Eigen::Vector3d const zz = v.cross(yy).normalized();
// y"_vec
- yy = MathLib::crossProduct(zz, v);
- yy.normalize();
+ Eigen::Vector3d const y = zz.cross(v).normalized();
- for (unsigned i=0; i<3; ++i) {
+ for (unsigned i = 0; i < 3; ++i)
+ {
rot_mat(0, i) = v[i];
- rot_mat(1, i) = yy[i];
+ rot_mat(1, i) = y[i];
rot_mat(2, i) = zz[i];
}
}
diff --git a/GeoLib/AnalyticalGeometry.cpp b/GeoLib/AnalyticalGeometry.cpp
index 5dc3bd94d5ca6a9012fab337bda3abe99eb3f5ed..cf481952d7ecfdf283dc9487673f3dcd8a52e904 100644
--- a/GeoLib/AnalyticalGeometry.cpp
+++ b/GeoLib/AnalyticalGeometry.cpp
@@ -377,17 +377,17 @@ std::unique_ptr<GeoLib::Point> triangleLineIntersection(
Eigen::Vector3d const pb = vb - vp;
Eigen::Vector3d const pc = vc - vp;
- double u(MathLib::scalarTriple(pq, pc, pb));
+ double u = pq.cross(pc).dot(pb);
if (u < 0)
{
return nullptr;
}
- double v(MathLib::scalarTriple(pq, pa, pc));
+ double v = pq.cross(pa).dot(pc);
if (v < 0)
{
return nullptr;
}
- double w(MathLib::scalarTriple(pq, pb, pa));
+ double w = pq.cross(pb).dot(pa);
if (w < 0)
{
return nullptr;
diff --git a/GeoLib/AnalyticalGeometry.h b/GeoLib/AnalyticalGeometry.h
index 5de5f0fa5cd1cf1beeb17f1619043a169c4ed938..34ce39f5e88b485311cc581d78896341cad18d34 100644
--- a/GeoLib/AnalyticalGeometry.h
+++ b/GeoLib/AnalyticalGeometry.h
@@ -86,7 +86,7 @@ std::pair<Eigen::Vector3d, double> getNewellPlane(
* @param rot_mat a 2x2 rotation matrix
*/
template<class T_MATRIX>
-void compute2DRotationMatrixToX(MathLib::Vector3 const& v, T_MATRIX & rot_mat);
+void compute2DRotationMatrixToX(Eigen::Vector3d const& v, T_MATRIX & rot_mat);
/**
* Computes a rotation matrix that rotates the given 3D normal vector parallel to X-axis
@@ -94,7 +94,7 @@ void compute2DRotationMatrixToX(MathLib::Vector3 const& v, T_MATRIX & rot_mat);
* @param rot_mat a 3x3 rotation matrix
*/
template <class T_MATRIX>
-void compute3DRotationMatrixToX(MathLib::Vector3 const& v, T_MATRIX& rot_mat);
+void compute3DRotationMatrixToX(Eigen::Vector3d const& v, T_MATRIX& rot_mat);
/**
* Method computes the rotation matrix that rotates the given vector parallel to
diff --git a/GeoLib/MinimalBoundingSphere.cpp b/GeoLib/MinimalBoundingSphere.cpp
index 115cc7facd20af4b668d8fac998322a486435f9b..64bac9e610f2a67888502bd13a7313495ae3f240 100644
--- a/GeoLib/MinimalBoundingSphere.cpp
+++ b/GeoLib/MinimalBoundingSphere.cpp
@@ -94,7 +94,7 @@ MinimalBoundingSphere::MinimalBoundingSphere(MathLib::Point3d const& p,
if (!MathLib::isCoplanar(p, q, r, s))
{
- double const denom = 2.0 * MathLib::scalarTriple(va, vb, vc);
+ double const denom = 2.0 * va.cross(vb).dot(vc);
Eigen::Vector3d o =
(vc.dot(vc) * va.cross(vb) + vb.dot(vb) * vc.cross(va) +
va.dot(va) * vb.cross(vc)) /
diff --git a/MathLib/GeometricBasics.cpp b/MathLib/GeometricBasics.cpp
index a9dbc23e341a251974e117b4b8da486c1728041a..0ce466a29bcfc87097def05d120db4072b64c946 100644
--- a/MathLib/GeometricBasics.cpp
+++ b/MathLib/GeometricBasics.cpp
@@ -12,7 +12,6 @@
#include <Eigen/Dense>
#include "Point3d.h"
-#include "Vector3.h"
#include "GeometricBasics.h"
@@ -28,10 +27,10 @@ double orientation3d(MathLib::Point3d const& p,
auto const pb = Eigen::Map<Eigen::Vector3d const>(b.getCoords());
auto const pc = Eigen::Map<Eigen::Vector3d const>(c.getCoords());
- Eigen::Vector3d const ap = pp - pa;
- Eigen::Vector3d const bp = pp - pb;
- Eigen::Vector3d const cp = pp - pc;
- return MathLib::scalarTriple(bp, cp, ap);
+ Eigen::Vector3d const u = pp - pa;
+ Eigen::Vector3d const v = pp - pb;
+ Eigen::Vector3d const w = pp - pc;
+ return u.cross(v).dot(w);
}
double calcTetrahedronVolume(MathLib::Point3d const& a,
@@ -43,19 +42,22 @@ double calcTetrahedronVolume(MathLib::Point3d const& a,
auto const vb = Eigen::Map<Eigen::Vector3d const>(b.getCoords());
auto const vc = Eigen::Map<Eigen::Vector3d const>(c.getCoords());
auto const vd = Eigen::Map<Eigen::Vector3d const>(d.getCoords());
- Eigen::Vector3d const ab = vb - va;
- Eigen::Vector3d const ac = vc - va;
- Eigen::Vector3d const ad = vd - va;
- return std::abs(MathLib::scalarTriple(ac, ad, ab)) / 6.0;
+ Eigen::Vector3d const w = vb - va;
+ Eigen::Vector3d const u = vc - va;
+ Eigen::Vector3d const v = vd - va;
+ return std::abs(u.cross(v).dot(w)) / 6.0;
}
double calcTriangleArea(MathLib::Point3d const& a, MathLib::Point3d const& b,
MathLib::Point3d 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();
+ auto const va = Eigen::Map<Eigen::Vector3d const>(a.getCoords());
+ auto const vb = Eigen::Map<Eigen::Vector3d const>(b.getCoords());
+ auto const vc = Eigen::Map<Eigen::Vector3d const>(c.getCoords());
+ Eigen::Vector3d const u = vc - va;
+ Eigen::Vector3d const v = vb - va;
+ Eigen::Vector3d const w = u.cross(v);
+ return 0.5 * w.norm();
}
bool isPointInTetrahedron(MathLib::Point3d const& p, MathLib::Point3d const& a,
@@ -123,17 +125,20 @@ bool gaussPointInTriangle(MathLib::Point3d const& q,
double eps_pnt_out_of_plane,
double eps_pnt_out_of_tri)
{
- MathLib::Vector3 const v(a, b);
- MathLib::Vector3 const w(a, c);
+ auto const pa = Eigen::Map<Eigen::Vector3d const>(a.getCoords());
+ auto const pb = Eigen::Map<Eigen::Vector3d const>(b.getCoords());
+ auto const pc = Eigen::Map<Eigen::Vector3d const>(c.getCoords());
+ Eigen::Vector3d const v = pb - pa;
+ Eigen::Vector3d const w = pc - pa;
Eigen::Matrix2d mat;
- mat(0, 0) = v.getSqrLength();
+ mat(0, 0) = v.squaredNorm();
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();
- Eigen::Vector2d y;
- y << 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]);
+ mat(1, 1) = w.squaredNorm();
+ Eigen::Vector2d y(
+ 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]));
y = mat.partialPivLu().solve(y);
@@ -167,17 +172,21 @@ bool barycentricPointInTriangle(MathLib::Point3d const& p,
return false;
}
- MathLib::Vector3 const pa(p, a);
- MathLib::Vector3 const pb(p, b);
- MathLib::Vector3 const pc(p, c);
+ auto const vp = Eigen::Map<Eigen::Vector3d const>(p.getCoords());
+ auto const va = Eigen::Map<Eigen::Vector3d const>(a.getCoords());
+ auto const vb = Eigen::Map<Eigen::Vector3d const>(b.getCoords());
+ auto const vc = Eigen::Map<Eigen::Vector3d const>(c.getCoords());
+ Eigen::Vector3d const pa = va - vp;
+ Eigen::Vector3d const pb = vb - vp;
+ Eigen::Vector3d const pc = vc - vp;
double const area_x_2(calcTriangleArea(a, b, c) * 2);
- double const alpha((MathLib::crossProduct(pb, pc).getLength()) / area_x_2);
+ double const alpha((pb.cross(pc).norm()) / area_x_2);
if (alpha < -eps_pnt_out_of_tri || alpha > 1 + eps_pnt_out_of_tri)
{
return false;
}
- double const beta((MathLib::crossProduct(pc, pa).getLength()) / area_x_2);
+ double const beta((pc.cross(pa).norm()) / area_x_2);
if (beta < -eps_pnt_out_of_tri || beta > 1 + eps_pnt_out_of_tri)
{
return false;
@@ -228,25 +237,30 @@ bool dividedByPlane(const MathLib::Point3d& a, const MathLib::Point3d& b,
bool isCoplanar(const MathLib::Point3d& a, const MathLib::Point3d& b,
const MathLib::Point3d& c, const MathLib::Point3d& d)
{
- const MathLib::Vector3 ab(a, b);
- const MathLib::Vector3 ac(a, c);
- const MathLib::Vector3 ad(a, d);
+ auto const pa = Eigen::Map<Eigen::Vector3d const>(a.getCoords());
+ auto const pb = Eigen::Map<Eigen::Vector3d const>(b.getCoords());
+ auto const pc = Eigen::Map<Eigen::Vector3d const>(c.getCoords());
+ auto const pd = Eigen::Map<Eigen::Vector3d const>(d.getCoords());
+
+ Eigen::Vector3d const ab = pb - pa;
+ Eigen::Vector3d const ac = pc - pa;
+ Eigen::Vector3d const ad = pd - pa;
- if (ab.getSqrLength() < pow(std::numeric_limits<double>::epsilon(), 2) ||
- ac.getSqrLength() < pow(std::numeric_limits<double>::epsilon(), 2) ||
- ad.getSqrLength() < pow(std::numeric_limits<double>::epsilon(), 2))
+ auto const eps_squared =
+ std::pow(std::numeric_limits<double>::epsilon(), 2);
+ if (ab.squaredNorm() < eps_squared || ac.squaredNorm() < eps_squared ||
+ ad.squaredNorm() < eps_squared)
{
return true;
}
// In exact arithmetic <ac*ad^T, ab> should be zero
// if all four points are coplanar.
- const double sqr_scalar_triple(
- pow(MathLib::scalarProduct(MathLib::crossProduct(ac, ad), ab), 2));
+ const double sqr_scalar_triple(std::pow(ac.cross(ad).dot(ab), 2));
// Due to evaluating the above numerically some cancellation or rounding
// can occur. For this reason a normalisation factor is introduced.
const double normalisation_factor =
- (ab.getSqrLength() * ac.getSqrLength() * ad.getSqrLength());
+ (ab.squaredNorm() * ac.squaredNorm() * ad.squaredNorm());
// tolerance 1e-11 is choosen such that
// a = (0,0,0), b=(1,0,0), c=(0,1,0) and d=(1,1,1e-6) are considered as
diff --git a/MathLib/MathTools.cpp b/MathLib/MathTools.cpp
index 5d7017408e0e27d0471fc812588c52c4f8b96509..ba12e31768e97b348630d1809863e21419062e45 100644
--- a/MathLib/MathTools.cpp
+++ b/MathLib/MathTools.cpp
@@ -49,14 +49,7 @@ double getAngle(Point3d const& p0, Point3d const& p1, Point3d const& p2)
Eigen::Vector3d const v1 = c - b;
// apply Cauchy Schwarz inequality
- return std::acos(
- (v0.transpose() * v1 / (v0.norm() * v1.norm()))(0, 0));
-}
-
-double scalarTriple(Eigen::Vector3d const& u, Eigen::Vector3d const& v,
- Eigen::Vector3d const& w)
-{
- return u.cross(v).dot(w);
+ return std::acos(v0.dot(v1) / (v0.norm() * v1.norm()));
}
} // namespace MathLib
diff --git a/MathLib/MathTools.h b/MathLib/MathTools.h
index 699d3ff7225ed3ea8c2368d25c29b90ad54060d2..812fb833fc922bb493594e4edd4392b05a58b889 100644
--- a/MathLib/MathTools.h
+++ b/MathLib/MathTools.h
@@ -91,8 +91,4 @@ double calcProjPntToLineAndDists(MathLib::Point3d const& pp,
* @return the angle between the edges
*/
double getAngle(Point3d const& p0, Point3d const& p1, Point3d const& p2);
-
-/// Calculates the scalar triple (u x v) . w
-double scalarTriple(Eigen::Vector3d const& u, Eigen::Vector3d const& v,
- Eigen::Vector3d const& w);
} // namespace MathLib
diff --git a/MeshGeoToolsLib/GeoMapper.cpp b/MeshGeoToolsLib/GeoMapper.cpp
index d3ed7fa7098b8b4207d3f353e300e7419751b1c6..41890c56780b293aeb16f33d9d080ec9b8c93952 100644
--- a/MeshGeoToolsLib/GeoMapper.cpp
+++ b/MeshGeoToolsLib/GeoMapper.cpp
@@ -33,7 +33,6 @@
#include "MeshLib/Elements/FaceRule.h"
#include "MeshLib/Node.h"
#include "MeshLib/MeshSurfaceExtraction.h"
-#include "MeshLib/MeshEditing/projectMeshOntoPlane.h"
#include "MeshLib/MeshSearch/MeshElementGrid.h"
namespace MeshGeoToolsLib {
@@ -419,12 +418,13 @@ static void mapPointOnSurfaceElement(MeshLib::Element const& elem,
MathLib::Point3d& q)
{
// create plane equation: n*p = d
- MathLib::Vector3 const p(*(elem.getNode(0)));
- MathLib::Vector3 const n(MeshLib::FaceRule::getSurfaceNormal(&elem));
+ auto const p =
+ Eigen::Map<Eigen::Vector3d const>(elem.getNode(0)->getCoords());
+ Eigen::Vector3d const n(MeshLib::FaceRule::getSurfaceNormal(&elem));
if (n[2] == 0.0) { // vertical plane, z coordinate is arbitrary
q[2] = p[2];
} else {
- double const d(MathLib::scalarProduct(n, p));
+ double const d(n.dot(p));
q[2] = (d - n[0]*q[0] - n[1]*q[1])/n[2];
}
}
diff --git a/MeshLib/CoordinateSystem.cpp b/MeshLib/CoordinateSystem.cpp
index ec26b6adfb57ef6deb446ced1db1d6fe4edaf481..83ace51cd3e83e0541254120285549289e88e340 100644
--- a/MeshLib/CoordinateSystem.cpp
+++ b/MeshLib/CoordinateSystem.cpp
@@ -41,7 +41,11 @@ unsigned char CoordinateSystem::getCoordinateSystem(const GeoLib::AABB &bbox) co
{
unsigned char coords = 0;
- const MathLib::Vector3 pt_diff(bbox.getMinPoint(), bbox.getMaxPoint());
+ auto const bbox_min =
+ Eigen::Map<Eigen::Vector3d const>(bbox.getMinPoint().getCoords());
+ auto const bbox_max =
+ Eigen::Map<Eigen::Vector3d const>(bbox.getMaxPoint().getCoords());
+ Eigen::Vector3d const pt_diff = bbox_max - bbox_min;
// The axis aligned bounding box is a from the right half-open interval.
// Therefore, the difference between the particular coordinates of the
diff --git a/MeshLib/CoordinateSystem.h b/MeshLib/CoordinateSystem.h
index 8bad010138a3c5c6b856bb1928607215578a4cd9..005111e26976f457b3f94abb70850cdf18bce1dc 100644
--- a/MeshLib/CoordinateSystem.h
+++ b/MeshLib/CoordinateSystem.h
@@ -12,7 +12,6 @@
#include <cmath>
#include "GeoLib/AABB.h"
-#include "MathLib/Vector3.h"
namespace MeshLib
{
diff --git a/MeshLib/ElementCoordinatesMappingLocal.cpp b/MeshLib/ElementCoordinatesMappingLocal.cpp
index 1d548a41d207fb5a4b1bd304c0e99b3cf46e630a..9b3500d9b1236dadaa87c1a03549c6fd60098f77 100644
--- a/MeshLib/ElementCoordinatesMappingLocal.cpp
+++ b/MeshLib/ElementCoordinatesMappingLocal.cpp
@@ -42,7 +42,9 @@ void getRotationMatrixToGlobal(const unsigned element_dimension,
// coordinates
if (element_dimension == 1)
{
- MathLib::Vector3 xx(points[0], points[1]);
+ auto const a = Eigen::Map<Eigen::Vector3d const>(points[0].getCoords());
+ auto const b = Eigen::Map<Eigen::Vector3d const>(points[1].getCoords());
+ Eigen::Vector3d xx = b - a;
xx.normalize();
if (global_dim == 2)
{
diff --git a/MeshLib/Elements/CellRule.cpp b/MeshLib/Elements/CellRule.cpp
index 5b10148112d0743f4bfdd728a0566051a8207147..3e5fe16fb9d5f70d06360e4f4112ec8e413fbef4 100644
--- a/MeshLib/Elements/CellRule.cpp
+++ b/MeshLib/Elements/CellRule.cpp
@@ -20,6 +20,7 @@ namespace MeshLib {
bool CellRule::testElementNodeOrder(const Element* e)
{
const MathLib::Vector3 c(getCenterOfGravity(*e));
+ Eigen::Vector3d const cc = Eigen::Map<Eigen::Vector3d const>(c.getCoords());
const unsigned nFaces (e->getNumberOfFaces());
for (unsigned j=0; j<nFaces; ++j)
{
@@ -28,9 +29,10 @@ bool CellRule::testElementNodeOrder(const Element* e)
// at some point, while for node 0 at least one node in every element
// type would be used for checking twice and one wouldn't be checked at
// all. (based on the definition of the _face_nodes variable)
- const MeshLib::Node x (*(face->getNode(1)));
- const MathLib::Vector3 cx (c, x);
- const double s = MathLib::scalarProduct(FaceRule::getSurfaceNormal(face), cx);
+ auto const x =
+ Eigen::Map<Eigen::Vector3d const>(face->getNode(1)->getCoords());
+ Eigen::Vector3d const cx = x - cc;
+ const double s = FaceRule::getSurfaceNormal(face).dot(cx);
delete face;
if (s >= 0)
{
diff --git a/MeshLib/Elements/FaceRule.cpp b/MeshLib/Elements/FaceRule.cpp
index 230c496db1cfdbccfdc405dc81fe13feb90c50b6..7d3d4e20e17ecf788ec1b3448163e299d8a82af5 100644
--- a/MeshLib/Elements/FaceRule.cpp
+++ b/MeshLib/Elements/FaceRule.cpp
@@ -21,23 +21,31 @@ bool FaceRule::testElementNodeOrder(const Element* e)
return getSurfaceNormal(e)[2] < 0;
}
-MathLib::Vector3 FaceRule::getFirstSurfaceVector(Element const* const e)
+Eigen::Vector3d FaceRule::getFirstSurfaceVector(Element const* const e)
{
- Node* const* const _nodes = e->getNodes();
- return {*_nodes[1], *_nodes[0]};
+ auto const a =
+ Eigen::Map<Eigen::Vector3d const>(e->getNode(0)->getCoords());
+ auto const b =
+ Eigen::Map<Eigen::Vector3d const>(e->getNode(1)->getCoords());
+ Eigen::Vector3d const v = a - b;
+ return v;
}
-MathLib::Vector3 FaceRule::getSecondSurfaceVector(Element const* const e)
+Eigen::Vector3d FaceRule::getSecondSurfaceVector(Element const* const e)
{
- Node* const* const _nodes = e->getNodes();
- return {*_nodes[1], *_nodes[2]};
+ auto const a =
+ Eigen::Map<Eigen::Vector3d const>(e->getNode(1)->getCoords());
+ auto const b =
+ Eigen::Map<Eigen::Vector3d const>(e->getNode(2)->getCoords());
+ Eigen::Vector3d const v = b - a;
+ return v;
}
-MathLib::Vector3 FaceRule::getSurfaceNormal(const Element* e)
+Eigen::Vector3d FaceRule::getSurfaceNormal(const Element* e)
{
- const MathLib::Vector3 u = getFirstSurfaceVector(e);
- const MathLib::Vector3 v = getSecondSurfaceVector(e);
- return MathLib::crossProduct(u, v);
+ Eigen::Vector3d const u = getFirstSurfaceVector(e);
+ Eigen::Vector3d const v = getSecondSurfaceVector(e);
+ return u.cross(v);
}
} // namespace MeshLib
diff --git a/MeshLib/Elements/FaceRule.h b/MeshLib/Elements/FaceRule.h
index e746a3d75fd3791e9adc3e1cae4100eab6c11c41..70d6e23e670eb13a24a5d304a4e8c9bd1a3b2112 100644
--- a/MeshLib/Elements/FaceRule.h
+++ b/MeshLib/Elements/FaceRule.h
@@ -35,13 +35,13 @@ public:
static bool testElementNodeOrder(const Element* /*e*/);
/// \returns the first vector forming the surface' plane
- static MathLib::Vector3 getFirstSurfaceVector(Element const* const e);
+ static Eigen::Vector3d getFirstSurfaceVector(Element const* const e);
/// \returns the second vector forming the surface' plane
- static MathLib::Vector3 getSecondSurfaceVector(Element const* const e);
+ static Eigen::Vector3d getSecondSurfaceVector(Element const* const e);
/// Returns the surface normal of a 2D element.
- static MathLib::Vector3 getSurfaceNormal(const Element* e);
+ static Eigen::Vector3d getSurfaceNormal(const Element* e);
}; /* class */
diff --git a/MeshLib/Elements/Utils.h b/MeshLib/Elements/Utils.h
index ae5c212c5223c4979dc85c46fd6a4e5162374cfa..3f992df76421dbcc099a48c88e0ebc59f76984b2 100644
--- a/MeshLib/Elements/Utils.h
+++ b/MeshLib/Elements/Utils.h
@@ -41,19 +41,21 @@ inline std::vector<Node*> getBaseNodes(std::vector<Element*> const& elements)
return base_nodes;
}
-inline MathLib::Vector3 calculateNormalizedSurfaceNormal(
+inline Eigen::Vector3d calculateNormalizedSurfaceNormal(
MeshLib::Element const& surface_element,
MeshLib::Element const& bulk_element)
{
- MathLib::Vector3 surface_element_normal;
+ Eigen::Vector3d surface_element_normal;
if (surface_element.getDimension() < 2)
{
- auto const bulk_element_normal = MeshLib::FaceRule::getSurfaceNormal(
- &bulk_element);
- MathLib::Vector3 const edge_vector(*surface_element.getNode(0),
- *surface_element.getNode(1));
- surface_element_normal =
- MathLib::crossProduct(bulk_element_normal, edge_vector);
+ auto const bulk_element_normal =
+ MeshLib::FaceRule::getSurfaceNormal(&bulk_element);
+ auto const v0 = Eigen::Map<Eigen::Vector3d const>(
+ surface_element.getNode(0)->getCoords());
+ auto const v1 = Eigen::Map<Eigen::Vector3d const>(
+ surface_element.getNode(1)->getCoords());
+ Eigen::Vector3d const edge_vector = v1 - v0;
+ surface_element_normal = bulk_element_normal.cross(edge_vector);
}
else
{
diff --git a/MeshLib/MeshEditing/ProjectPointOnMesh.cpp b/MeshLib/MeshEditing/ProjectPointOnMesh.cpp
index 0453a133d432bf5afb76984faf51614b2b03d29c..9a4e225f5a6ad56788ee06d930f4abca24fb14c6 100644
--- a/MeshLib/MeshEditing/ProjectPointOnMesh.cpp
+++ b/MeshLib/MeshEditing/ProjectPointOnMesh.cpp
@@ -58,14 +58,13 @@ MeshLib::Element const* getProjectedElement(
return nullptr;
}
-double getElevation(MeshLib::Element const& element,
- MeshLib::Node const& node)
+double getElevation(MeshLib::Element const& element, MeshLib::Node const& node)
{
- MathLib::Vector3 const v =
- MathLib::Vector3(node) - MathLib::Vector3(*element.getNode(0));
- MathLib::Vector3 const n =
- MeshLib::FaceRule::getSurfaceNormal(&element).getNormalizedVector();
- return node[2] - scalarProduct(n, v) * n[2];
+ Eigen::Vector3d const v =
+ Eigen::Map<Eigen::Vector3d const>(node.getCoords()) -
+ Eigen::Map<Eigen::Vector3d const>(element.getNode(0)->getCoords());
+ auto const n = MeshLib::FaceRule::getSurfaceNormal(&element).normalized();
+ return node[2] - n.dot(v) * n[2];
}
} // namespace ProjectPointOnMesh
diff --git a/MeshLib/MeshEditing/projectMeshOntoPlane.h b/MeshLib/MeshEditing/projectMeshOntoPlane.h
deleted file mode 100644
index ae269bf49572c734591f547869526e4be001bf79..0000000000000000000000000000000000000000
--- a/MeshLib/MeshEditing/projectMeshOntoPlane.h
+++ /dev/null
@@ -1,58 +0,0 @@
-/**
- * \file
- * \author Karsten Rink
- * \date 2015-04-10
- * \brief Definition of the projectMeshOntoPlane
- *
- * \copyright
- * Copyright (c) 2012-2020, OpenGeoSys Community (http://www.opengeosys.org)
- * Distributed under a Modified BSD License.
- * See accompanying file LICENSE.txt or
- * http://www.opengeosys.org/project/license
- *
- */
-
-#pragma once
-
-#include <vector>
-
-#include "MathLib/MathTools.h"
-#include "MathLib/Point3d.h"
-#include "MathLib/Vector3.h"
-
-#include "MeshLib/Mesh.h"
-#include "MeshLib/Node.h"
-#include "MeshLib/MeshEditing/DuplicateMeshComponents.h"
-
-
-namespace MeshLib {
-
-/**
- * Projects all nodes of a mesh onto a plane specified by a point of origin and a normal vector.
- * Overlapping elements, collapsed nodes, and other issues are not handled by the method.
- * The normal vector need not be normalized.
- */
-inline
-MeshLib::Mesh* projectMeshOntoPlane(MeshLib::Mesh const& mesh,
- MathLib::Point3d const& plane_origin,
- MathLib::Vector3 const& plane_normal)
-{
- std::size_t const n_nodes (mesh.getNumberOfNodes());
- std::vector<MeshLib::Node*> const& nodes (mesh.getNodes());
- MathLib::Vector3 normal (plane_normal);
- normal.normalize();
- std::vector<MeshLib::Node*> new_nodes;
- new_nodes.reserve(n_nodes);
- for (std::size_t i=0; i<n_nodes; ++i)
- {
- MeshLib::Node const& node(*nodes[i]);
- MathLib::Vector3 const v(plane_origin, node);
- double const dist (MathLib::scalarProduct(v,normal));
- new_nodes.push_back(new MeshLib::Node(node - dist * normal));
- }
-
- return new MeshLib::Mesh("Projected_Mesh", new_nodes,
- MeshLib::copyElementVector(mesh.getElements(), new_nodes));
-}
-
-} // end namespace MeshLib
diff --git a/MeshLib/MeshSurfaceExtraction.cpp b/MeshLib/MeshSurfaceExtraction.cpp
index 0f6996e78d6761078f6bb5f89c2ecc337a8dc2fd..76f792587ac5a29d6684693349fb8cab99a78901 100644
--- a/MeshLib/MeshSurfaceExtraction.cpp
+++ b/MeshLib/MeshSurfaceExtraction.cpp
@@ -289,6 +289,8 @@ void MeshSurfaceExtraction::get2DSurfaceElements(
std::vector<std::size_t>& element_to_bulk_face_id_map,
const MathLib::Vector3& dir, double angle, unsigned mesh_dimension)
{
+ auto const d = Eigen::Map<Eigen::Vector3d const>(dir.getCoords());
+
if (mesh_dimension < 2 || mesh_dimension > 3)
{
ERR("Cannot handle meshes of dimension {:i}", mesh_dimension);
@@ -298,7 +300,7 @@ void MeshSurfaceExtraction::get2DSurfaceElements(
double const pi(boost::math::constants::pi<double>());
double const cos_theta(std::cos(angle * pi / 180.0));
- MathLib::Vector3 const norm_dir(dir.getNormalizedVector());
+ Eigen::Vector3d const norm_dir(d.normalized());
for (auto const* elem : all_elements)
{
@@ -313,8 +315,7 @@ void MeshSurfaceExtraction::get2DSurfaceElements(
if (!complete_surface)
{
auto const* face = elem;
- if (MathLib::scalarProduct(
- FaceRule::getSurfaceNormal(face).getNormalizedVector(),
+ if (FaceRule::getSurfaceNormal(face).normalized().dot(
norm_dir) > cos_theta)
{
continue;
@@ -342,10 +343,9 @@ void MeshSurfaceExtraction::get2DSurfaceElements(
std::unique_ptr<MeshLib::Element const>{elem->getFace(j)};
if (!complete_surface)
{
- if (MathLib::scalarProduct(
- FaceRule::getSurfaceNormal(face.get())
- .getNormalizedVector(),
- norm_dir) < cos_theta)
+ if (FaceRule::getSurfaceNormal(face.get())
+ .normalized()
+ .dot(norm_dir) < cos_theta)
{
continue;
}
diff --git a/MeshLib/Node.h b/MeshLib/Node.h
index 2ee62b82e039a256d3f7df836a636fa1423a6e09..4a3f465321c5cbe4f690f47f497a897c58e826d3 100644
--- a/MeshLib/Node.h
+++ b/MeshLib/Node.h
@@ -69,12 +69,6 @@ public:
/// Get number of elements the node is part of.
std::size_t getNumberOfElements() const { return _elements.size(); }
- /// Shift the node according to the displacement vector v.
- Node operator-(MathLib::Vector3 const& v) const
- {
- return Node(_x[0]-v[0], _x[1]-v[1], _x[2]-v[2]);
- }
-
protected:
/// Update coordinates of a node.
/// This method automatically also updates the areas/volumes of all connected elements.
diff --git a/ProcessLib/BoundaryCondition/ConstraintDirichletBoundaryConditionLocalAssembler.h b/ProcessLib/BoundaryCondition/ConstraintDirichletBoundaryConditionLocalAssembler.h
index e456bab90a0344c92f21c89f538951dec354c162..af6608fc03dd158c981665d9f7ea09f2d169b1cf 100644
--- a/ProcessLib/BoundaryCondition/ConstraintDirichletBoundaryConditionLocalAssembler.h
+++ b/ProcessLib/BoundaryCondition/ConstraintDirichletBoundaryConditionLocalAssembler.h
@@ -132,8 +132,7 @@ public:
double const bulk_grad_times_normal(
Eigen::Map<Eigen::RowVectorXd const>(bulk_flux.data(),
bulk_flux.size())
- .dot(Eigen::Map<Eigen::RowVectorXd const>(
- _surface_element_normal.getCoords(), 3)));
+ .dot(_surface_element_normal));
integrated_value +=
bulk_grad_times_normal *
@@ -149,7 +148,7 @@ private:
IntegrationMethod const _integration_method;
std::size_t const _bulk_element_id;
- MathLib::Vector3 const _surface_element_normal;
+ Eigen::Vector3d const _surface_element_normal;
};
} // namespace ProcessLib
diff --git a/ProcessLib/BoundaryCondition/HCNonAdvectiveFreeComponentFlowBoundaryConditionLocalAssembler.h b/ProcessLib/BoundaryCondition/HCNonAdvectiveFreeComponentFlowBoundaryConditionLocalAssembler.h
index 9a67ffdbac4d3879264f3649cc5da9de9bc00860..e0e91258c2d455152b4fb0fa9baea41d61f9bbcb 100644
--- a/ProcessLib/BoundaryCondition/HCNonAdvectiveFreeComponentFlowBoundaryConditionLocalAssembler.h
+++ b/ProcessLib/BoundaryCondition/HCNonAdvectiveFreeComponentFlowBoundaryConditionLocalAssembler.h
@@ -54,18 +54,6 @@ public:
_surface_normal(getOrientedSurfaceNormal(e))
{
}
- Eigen::RowVector3d getOrientedSurfaceNormal(MeshLib::Element const& e)
- {
- auto surface_element_normal = MeshLib::FaceRule::getSurfaceNormal(&e);
- surface_element_normal.normalize();
- // At the moment (2016-09-28) the surface normal is not oriented
- // according to the right hand rule
- // for correct results it is necessary to multiply the normal with -1
- surface_element_normal *= -1;
- return Eigen::Map<Eigen::RowVector3d const>(
- surface_element_normal.getCoords(),
- _data.process.getMesh().getDimension());
- }
void assemble(std::size_t const mesh_item_id,
NumLib::LocalToGlobalIndexMap const& dof_table_boundary,
@@ -113,9 +101,21 @@ public:
}
private:
+ Eigen::Vector3d getOrientedSurfaceNormal(MeshLib::Element const& e)
+ {
+ // At the moment (2016-09-28) the surface normal is not oriented
+ // according to the right hand rule
+ // for correct results it is necessary to multiply the normal with -1
+ Eigen::Vector3d surface_normal =
+ -MeshLib::FaceRule::getSurfaceNormal(&e).normalized();
+ auto const zeros_size = 3 - _data.process.getMesh().getDimension();
+ surface_normal.tail(zeros_size).setZero();
+ return surface_normal;
+ }
+
HCNonAdvectiveFreeComponentFlowBoundaryConditionData const& _data;
std::size_t const _local_matrix_size;
- Eigen::RowVector3d const _surface_normal;
+ Eigen::Vector3d const _surface_normal;
};
} // namespace ProcessLib
diff --git a/ProcessLib/BoundaryCondition/NormalTractionBoundaryConditionLocalAssembler.h b/ProcessLib/BoundaryCondition/NormalTractionBoundaryConditionLocalAssembler.h
index d32103068b7f7a8bf82dc700b64cee0f52bfb104..47f5b685519297130c4ac7f8cd64a7959e9598ad 100644
--- a/ProcessLib/BoundaryCondition/NormalTractionBoundaryConditionLocalAssembler.h
+++ b/ProcessLib/BoundaryCondition/NormalTractionBoundaryConditionLocalAssembler.h
@@ -94,11 +94,11 @@ public:
}
else
{
- auto const element_normal_vector =
- MeshLib::FaceRule::getSurfaceNormal(&e).getNormalizedVector();
-
- std::copy_n(element_normal_vector.getCoords(), GlobalDim,
- element_normal.data());
+ auto const n = MeshLib::FaceRule::getSurfaceNormal(&e).normalized();
+ for (decltype(GlobalDim) i = 0; i < GlobalDim; ++i)
+ {
+ element_normal[i] = n[i];
+ }
}
for (unsigned ip = 0; ip < n_integration_points; ip++)
diff --git a/ProcessLib/LIE/Common/Utils.cpp b/ProcessLib/LIE/Common/Utils.cpp
index d6c52f604c17737b6f1ed29933118ff9fc7d6053..f18601f1c78b2373deb5b08ebae2501652ca6223 100644
--- a/ProcessLib/LIE/Common/Utils.cpp
+++ b/ProcessLib/LIE/Common/Utils.cpp
@@ -15,6 +15,7 @@ namespace ProcessLib
{
namespace LIE
{
+// ToDo (TF) change interface
void computeNormalVector(MeshLib::Element const& e, unsigned const global_dim,
Eigen::Vector3d& element_normal)
{
@@ -31,10 +32,11 @@ void computeNormalVector(MeshLib::Element const& e, unsigned const global_dim,
else if (global_dim == 3)
{
auto const element_normal_vector =
- MeshLib::FaceRule::getSurfaceNormal(&e).getNormalizedVector();
-
- std::copy_n(element_normal_vector.getCoords(), global_dim,
- element_normal.data());
+ MeshLib::FaceRule::getSurfaceNormal(&e).normalized();
+ for (int i = 0; i < 3; ++i)
+ {
+ element_normal[i] = element_normal_vector[i];
+ }
}
}
@@ -49,9 +51,9 @@ void computeRotationMatrix(MeshLib::Element const& e, Eigen::Vector3d const& n,
else if (global_dim == 3)
{
auto const u =
- MeshLib::FaceRule::getFirstSurfaceVector(&e).getNormalizedVector();
+ MeshLib::FaceRule::getFirstSurfaceVector(&e).normalized();
auto const v =
- MeshLib::FaceRule::getSecondSurfaceVector(&e).getNormalizedVector();
+ MeshLib::FaceRule::getSecondSurfaceVector(&e).normalized();
R.resize(3, 3);
R << u[0], u[1], u[2], v[0], v[1], v[2], n[0], n[1], n[2];
diff --git a/ProcessLib/SurfaceFlux/SurfaceFluxLocalAssembler.h b/ProcessLib/SurfaceFlux/SurfaceFluxLocalAssembler.h
index f1de6d115c13ee331e3ec533a7b8bac3d633a007..82654e6d10bc7a1fda52134c90a23fe7f4cdc6c1 100644
--- a/ProcessLib/SurfaceFlux/SurfaceFluxLocalAssembler.h
+++ b/ProcessLib/SurfaceFlux/SurfaceFluxLocalAssembler.h
@@ -109,16 +109,18 @@ public:
{
auto get_surface_normal =
[this, &bulk_mesh](
- MeshLib::Element const& surface_element) -> MathLib::Vector3 {
- MathLib::Vector3 surface_element_normal;
+ MeshLib::Element const& surface_element) -> Eigen::Vector3d {
+ Eigen::Vector3d surface_element_normal;
if (surface_element.getGeomType() == MeshLib::MeshElemType::LINE)
{
auto const bulk_normal = MeshLib::FaceRule::getSurfaceNormal(
bulk_mesh.getElements()[_bulk_element_id]);
- MathLib::Vector3 const line{*_surface_element.getNodes()[0],
- *_surface_element.getNodes()[1]};
- surface_element_normal =
- MathLib::crossProduct(line, bulk_normal);
+ auto const l0 = Eigen::Map<Eigen::Vector3d const>(
+ _surface_element.getNode(0)->getCoords());
+ auto const l1 = Eigen::Map<Eigen::Vector3d const>(
+ _surface_element.getNode(1)->getCoords());
+ Eigen::Vector3d const line = l1 - l0;
+ surface_element_normal = line.cross(bulk_normal);
}
else
{
@@ -157,8 +159,7 @@ public:
double const bulk_grad_times_normal(
Eigen::Map<Eigen::RowVectorXd const>(bulk_flux.data(),
bulk_flux.size())
- .dot(Eigen::Map<Eigen::RowVectorXd const>(
- surface_element_normal.getCoords(), 3)));
+ .dot(surface_element_normal));
specific_flux.getComponent(element_id, component_id) +=
bulk_grad_times_normal * _detJ_times_integralMeasure[ip] *
diff --git a/Tests/MeshLib/TestProjectMeshOnPlane.cpp b/Tests/MeshLib/TestProjectMeshOnPlane.cpp
deleted file mode 100644
index 4cc080f1efa1de8437b4b0f245d1406c9c44e3ea..0000000000000000000000000000000000000000
--- a/Tests/MeshLib/TestProjectMeshOnPlane.cpp
+++ /dev/null
@@ -1,140 +0,0 @@
-/**
- * \file
- * \author Karsten Rink
- * \date 2015-04-16
- * \brief Tests for projectMeshOnPlane
- *
- * \copyright
- * Copyright (c) 2012-2020, 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 <memory>
-
-#include "gtest/gtest.h"
-
-#include "MeshLib/Elements/Element.h"
-#include "MeshLib/Elements/Line.h"
-#include "MeshLib/Mesh.h"
-#include "MeshLib/MeshEditing/projectMeshOntoPlane.h"
-#include "MeshLib/Node.h"
-
-class ProjectionTest : public ::testing::Test
-{
-public:
- ProjectionTest()
- : _mesh(nullptr)
- {
- std::size_t const n_nodes (100);
- std::vector<MeshLib::Node*> nodes;
- for (std::size_t i = 1; i <= n_nodes; i++)
- {
- nodes.push_back(new MeshLib::Node(static_cast<double>(i),
- static_cast<double>(i), i * 0.5));
- }
-
- std::vector<MeshLib::Element*> elements;
- for (std::size_t i = 0; i < n_nodes - 1; i++)
- {
- elements.push_back(new MeshLib::Line(
- std::array<MeshLib::Node*, 2>{{nodes[i], nodes[i + 1]}}));
- }
-
- _mesh = std::make_unique<MeshLib::Mesh>("TestMesh", nodes, elements);
- }
-
-protected:
- std::unique_ptr<MeshLib::Mesh> _mesh;
-};
-// Project to parallels of XY plane
-TEST_F(ProjectionTest, ProjectToXY)
-{
- MathLib::Vector3 normal (0,0,1);
- std::size_t const n_nodes (_mesh->getNumberOfNodes());
- for (std::size_t p=0; p<10; p++)
- {
- MathLib::Point3d origin (std::array<double,3>{{0,0,static_cast<double>(p)}});
- MeshLib::Mesh* result = MeshLib::projectMeshOntoPlane(*_mesh, origin, normal);
- for (std::size_t i = 0; i < n_nodes; i++)
- {
- ASSERT_NEAR(static_cast<double>(p), (*result->getNode(i))[2],
- std::numeric_limits<double>::epsilon());
- }
- delete result;
- }
-}
-
-// Project to parallels of XZ plane
-TEST_F(ProjectionTest, ProjectToXZ)
-{
- MathLib::Vector3 normal (0,1,0);
- std::size_t const n_nodes (_mesh->getNumberOfNodes());
- for (std::size_t p=0; p<10; p++)
- {
- MathLib::Point3d origin (std::array<double,3>{{0,static_cast<double>(p),0}});
- MeshLib::Mesh* result = MeshLib::projectMeshOntoPlane(*_mesh, origin, normal);
- for (std::size_t i = 0; i < n_nodes; i++)
- {
- ASSERT_NEAR(static_cast<double>(p), (*result->getNode(i))[1],
- std::numeric_limits<double>::epsilon());
- }
- delete result;
- }
-}
-
-// Project to parallels of YZ plane
-TEST_F(ProjectionTest, ProjectToYZ)
-{
- MathLib::Vector3 normal (1,0,0);
- std::size_t const n_nodes (_mesh->getNumberOfNodes());
- for (std::size_t p=0; p<10; p++)
- {
- MathLib::Point3d origin (std::array<double,3>{{static_cast<double>(p),0,0}});
- MeshLib::Mesh* result = MeshLib::projectMeshOntoPlane(*_mesh, origin, normal);
- for (std::size_t i = 0; i < n_nodes; i++)
- {
- ASSERT_NEAR(static_cast<double>(p), (*result->getNode(i))[0],
- std::numeric_limits<double>::epsilon());
- }
- delete result;
- }
-}
-
-// Sign of normal vector does not matter.
-TEST_F(ProjectionTest, NormalDirection)
-{
- MathLib::Vector3 normal_p (0,0,1);
- MathLib::Vector3 normal_n (0,0,-1);
- std::size_t const n_nodes (_mesh->getNumberOfNodes());
- MathLib::Point3d origin (std::array<double,3>{{0,0,0}});
- MeshLib::Mesh* result_p = MeshLib::projectMeshOntoPlane(*_mesh, origin, normal_p);
- MeshLib::Mesh* result_n = MeshLib::projectMeshOntoPlane(*_mesh, origin, normal_n);
- for (std::size_t i = 0; i < n_nodes; i++)
- {
- ASSERT_EQ((*result_p->getNode(i))[2], (*result_n->getNode(i))[2]);
- }
- delete result_p;
- delete result_n;
-}
-
-// Length of normal does not matter (it's normalised within the method)
-TEST_F(ProjectionTest, NormalLength)
-{
- MathLib::Point3d origin (std::array<double,3>{{0,0,0}});
- MathLib::Vector3 normal (0,0,1);
- std::size_t const n_nodes (_mesh->getNumberOfNodes());
- MeshLib::Mesh* result = MeshLib::projectMeshOntoPlane(*_mesh, origin, normal);
- for (std::size_t p=2; p<10; p++)
- {
- normal[2] = static_cast<double>(p);
- MeshLib::Mesh* result_p = MeshLib::projectMeshOntoPlane(*_mesh, origin, normal);
- for (std::size_t i = 0; i < n_nodes; i++)
- {
- ASSERT_EQ((*result->getNode(i))[2], (*result_p->getNode(i))[2]);
- }
- delete result_p;
- }
- delete result;
-}