OGS
Go to the documentation of this file.
1
10
11
12#include <Eigen/Geometry>
13#include <algorithm>
14#include <memory>
15#include <string>
16#include <vector>
17
18#include "GeoLib/AABB.h"
19#include "InfoLib/GitInfo.h"
20#include "MathLib/Point3d.h"
24#include "MeshLib/Mesh.h"
25#include "MeshLib/Node.h"
26
27namespace
28{
29// tests if a plane and an AABB are intersecting
30// (based on Christer Ericson "Real Time Collision Detection" 5.2.3)
31bool testAABBIntersectingPlane(Eigen::Vector3d const& aabb_centre,
32 Eigen::Vector3d const& aabb_extent,
33 Eigen::Vector3d const& plane_normal,
34 double const pd)
35{
36 double const r = aabb_extent.dot(plane_normal.cwiseAbs());
37 double const s = plane_normal.dot(aabb_centre) - pd;
38 return std::abs(s) > r;
39}
40// tests if a triangle and an AABB are intersecting
41// (based on Christer Ericson "Real Time Collision Detection" 5.2.9)
43 MeshLib::Node const& n1,
44 MeshLib::Node const& n2,
45 Eigen::Vector3d const& c,
46 Eigen::Vector3d const& e)
47{
48 // Translate triangle as conceptually moving AABB to origin
49 Eigen::Matrix3d v;
50 v << n0.asEigenVector3d() - c, n1.asEigenVector3d() - c,
51 n2.asEigenVector3d() - c;
52
53 // Test the three axes corresponding to the face normals of AABB b
54 if (((v.rowwise().minCoeff() - e).array() > 0).any() ||
55 ((v.rowwise().maxCoeff() + e).array() < 0).any())
56 {
57 return false;
58 }
59
60 // separating axes
61 std::array<Eigen::Vector3d, 3> tri_edge{
62 {v.col(1) - v.col(0), v.col(2) - v.col(1), v.col(0) - v.col(2)}};
63 std::array<Eigen::Vector3d, 9> const axx{
64 {Eigen::Vector3d({0, -tri_edge[0].z(), tri_edge[0].y()}),
65 Eigen::Vector3d({0, -tri_edge[1].z(), tri_edge[1].y()}),
66 Eigen::Vector3d({0, -tri_edge[2].z(), tri_edge[2].y()}),
67 Eigen::Vector3d({tri_edge[0].z(), 0, -tri_edge[0].x()}),
68 Eigen::Vector3d({tri_edge[1].z(), 0, -tri_edge[1].x()}),
69 Eigen::Vector3d({tri_edge[2].z(), 0, -tri_edge[2].x()}),
70 Eigen::Vector3d({-tri_edge[0].y(), tri_edge[0].x(), 0}),
71 Eigen::Vector3d({-tri_edge[1].y(), tri_edge[1].x(), 0}),
72 Eigen::Vector3d({-tri_edge[2].y(), tri_edge[2].x(), 0})}};
73
74 // Separating axis tests to check if there's a plane separating the
75 // projections of the AABB and the triangle according to the Separating Axis
76 // Theorem (see C. Ericson "Real Time Collision Detection" for details)
77 for (auto const& a : axx)
78 {
79 Eigen::Vector3d p = v.transpose() * a;
80 double const r = e.dot(a.cwiseAbs());
81 if (std::max(-p.maxCoeff(), p.minCoeff()) > r)
82 {
83 return false;
84 }
85 }
86
87 // Test separating axis corresponding to triangle face normal
88 Eigen::Vector3d const plane_normal(tri_edge[0].cross(tri_edge[1]));
89 double const pd = plane_normal.dot(v.row(0));
90 return testAABBIntersectingPlane(c, e, plane_normal, pd);
91}
92// mark all cells of the voxel grid that intersect with fault
93void markFaults(MeshLib::Mesh& mesh, MeshLib::Mesh const& fault,
94 int const fault_id, Eigen::Vector3d const& half_cell_size)
95{
96 auto const& elems = mesh.getElements();
97 std::size_t const n_elems = mesh.getNumberOfElements();
98 auto mat_ids = MeshLib::materialIDs(mesh);
99 auto const& fnodes = fault.getNodes();
100 auto const& felems = fault.getElements();
101 GeoLib::AABB const fault_aabb(fnodes.cbegin(), fnodes.cend());
102 auto [min_pnt, max_pnt] = fault_aabb.getMinMaxPoints();
103
104 // get bounding box of fault + voxel extent
105 min_pnt -= half_cell_size;
106 max_pnt += half_cell_size;
107
108 std::array<Eigen::Vector3d, 2> const fault_extent{{min_pnt, max_pnt}};
109 GeoLib::AABB const fault_aabb_ext(fault_extent.cbegin(),
110 fault_extent.cend());
111
112 // test each voxel grid element vs each fault triangle
113 Eigen::Vector3d const extent{half_cell_size};
114 for (std::size_t j = 0; j < n_elems; ++j)
115 {
116 // test if bounding box of fault is intersecting voxel
117 auto const& centre_pnt = MeshLib::getCenterOfGravity(*elems[j]);
118 if (!fault_aabb_ext.containsPoint(centre_pnt, 0))
119 {
120 continue;
121 }
122
123 // test if voxel is intersecting triangle
124 auto const& c(centre_pnt.asEigenVector3d());
125 for (auto const* const fault_elem : felems)
126 {
127 if (fault_elem->getDimension() != 2)
128 {
129 continue;
130 }
131
133 *fault_elem->getNode(0), *fault_elem->getNode(1),
134 *fault_elem->getNode(2), c, extent))
135 {
136 (*mat_ids)[j] = fault_id;
137 break;
138 }
139
141 {
143 *fault_elem->getNode(0), *fault_elem->getNode(2),
144 *fault_elem->getNode(3), c, extent))
145 {
146 (*mat_ids)[j] = fault_id;
147 break;
148 }
149 }
150 }
151 }
152}
153} // namespace
155{
156bool isVoxelGrid(MeshLib::Mesh const& mesh)
157{
158 auto const& elements = mesh.getElements();
159 if (std::any_of(elements.cbegin(), elements.cend(),
160 [&](auto const& e) {
161 return (e->getGeomType() !=
162 MeshLib::MeshElemType::HEXAHEDRON);
163 }))
164 {
165 return false;
166 }
167 auto is_voxel = [](auto const& e)
168 {
169 auto const n = e->getNodes();
170 return ((*n[0])[2] != (*n[1])[2] || (*n[1])[2] != (*n[2])[2] ||
171 (*n[4])[2] != (*n[5])[2] || (*n[5])[2] != (*n[6])[2] ||
172 (*n[1])[0] != (*n[2])[0] || (*n[2])[0] != (*n[5])[0] ||
173 (*n[0])[0] != (*n[3])[0] || (*n[3])[0] != (*n[7])[0]);
174 };
175
176 if (std::any_of(elements.cbegin(), elements.cend(), is_voxel))
177 {
178 ERR("Input mesh needs to be voxel grid (i.e. equally sized axis "
179 "aligned hexahedra).");
180 return false;
181 }
182 return true;
183}
186{
188 MeshLib::Mesh const* fault,
189 int const fault_id)
190{
191 if (mesh == nullptr)
192 {
194 return false;
195 }
196 if (!isVoxelGrid(*mesh))
197 {
198 ERR("The input mesh is not a voxel grid. The input mesh must be "
199 "a voxel grid (i.e. an equally sized axis "
200 "aligned hexahedra mesh).");
201 return false;
202 }
203
204 if (fault == nullptr)
205 {
207 return false;
208 }
209 if (fault->getDimension() != 2)
210 {
211 ERR("Fault needs to be a 2D mesh.");
212 return false;
213 }
214
215 Eigen::Vector3d half_cell_size;
216 {
217 auto const n = *mesh->getElement(0)->getNode(0);
218 auto const c = MeshLib::getCenterOfGravity(*mesh->getElement(0));
219 half_cell_size[0] = std::abs(c[0] - n[0]);
220 half_cell_size[1] = std::abs(c[1] - n[1]);
221 half_cell_size[2] = std::abs(c[2] - n[2]);
222 }
223
224 markFaults(*mesh, *fault, fault_id, half_cell_size);
225
226 return true;
227}
Definition of the AABB class.
Definition of the Element class.
Git information.
void ERR(fmt::format_string< Args... > fmt, Args &&... args)
Definition Logging.h:45
Definition of the Mesh class.
Definition of the Node class.
Definition of the Point3d class.
Implementation of the VtuInterface class.
Class AABB is an axis aligned bounding box around a given set of geometric points of (template) type ...
Definition AABB.h:56
bool containsPoint(T const &pnt, double eps) const
Definition AABB.h:143
MinMaxPoints getMinMaxPoints() const
Definition AABB.h:174
Eigen::Vector3d const & asEigenVector3d() const
Definition Point3d.h:64
virtual const Node * getNode(unsigned idx) const =0
std::vector< Node * > const & getNodes() const
Get the nodes-vector for the mesh.
Definition Mesh.h:106
std::vector< Element * > const & getElements() const
Get the element-vector for the mesh.
Definition Mesh.h:109
unsigned getDimension() const
Returns the dimension of the mesh (determined by the maximum dimension over all elements).
Definition Mesh.h:88
const Element * getElement(std::size_t idx) const
Get the element with the given index.
Definition Mesh.h:94
std::size_t getNumberOfElements() const
Get the number of elements.
Definition Mesh.h:97
PropertyVector< int > const * materialIDs(Mesh const &mesh)
Definition Mesh.cpp:268
MathLib::Point3d getCenterOfGravity(Element const &element)
Calculates the center of gravity for the mesh element.
Definition Element.cpp:124
bool addFaultToVoxelGrid(MeshLib::Mesh *mesh, MeshLib::Mesh const *fault, int const fault_id)
bool testAABBIntersectingPlane(Eigen::Vector3d const &aabb_centre, Eigen::Vector3d const &aabb_extent, Eigen::Vector3d const &plane_normal, double const pd)
void markFaults(MeshLib::Mesh &mesh, MeshLib::Mesh const &fault, int const fault_id, Eigen::Vector3d const &half_cell_size)
bool testTriangleIntersectingAABB(MeshLib::Node const &n0, MeshLib::Node const &n1, MeshLib::Node const &n2, Eigen::Vector3d const &c, Eigen::Vector3d const &e)