Skip to content

Commit

Permalink
Merge branch 'VoxelFault' into 'master'
Browse files Browse the repository at this point in the history 
[tools] Extract functions from AddFaultToVoxelGrid

See merge request ogs/ogs!4681
  • Loading branch information
@endJunction
endJunction committed Jul 26, 2023
2 parents 43a4b24 + e67fe83 commit 86515bc
Show file tree
Hide file tree
Showing 5 changed files with 275 additions and 206 deletions.
222 changes: 16 additions & 206 deletions Applications/Utils/MeshEdit/AddFaultToVoxelGrid.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -15,179 +15,15 @@
// ThirdParty
#include <tclap/CmdLine.h>

#include <Eigen/Geometry>

#ifdef USE_PETSC
#include <mpi.h>
#endif

#include "GeoLib/AABB.h"
#include "InfoLib/GitInfo.h"
#include "MathLib/Point3d.h"
#include "MeshLib/Elements/Element.h"
#include "MeshLib/IO/VtkIO/VtuInterface.h"
#include "MeshLib/IO/readMeshFromFile.h"
#include "MeshLib/Mesh.h"
#include "MeshLib/Node.h"

static std::string mat_name = "MaterialIDs";

// tests if a plane and an AABB are intersecting
// (based on Christer Ericson "Real Time Collision Detection" 5.2.3)
bool testAABBIntersectingPlane(Eigen::Vector3d const& aabb_centre,
Eigen::Vector3d const& aabb_extent,
Eigen::Vector3d const& plane_normal,
double const pd)
{
double const r = aabb_extent.dot(plane_normal.cwiseAbs());
double const s = plane_normal.dot(aabb_centre) - pd;
return std::abs(s) > r;
}

// tests if a triangle and an AABB are intersecting
// (based on Christer Ericson "Real Time Collision Detection" 5.2.9)
bool testTriangleIntersectingAABB(MeshLib::Node const& n0,
MeshLib::Node const& n1,
MeshLib::Node const& n2,
Eigen::Vector3d const& c,
Eigen::Vector3d const& e)
{
// Translate triangle as conceptually moving AABB to origin
Eigen::Matrix3d v;
v << n0.asEigenVector3d() - c, n1.asEigenVector3d() - c,
n2.asEigenVector3d() - c;

// Test the three axes corresponding to the face normals of AABB b
if (((v.rowwise().minCoeff() - e).array() > 0).any() ||
((v.rowwise().maxCoeff() + e).array() < 0).any())
{
return false;
}

// separating axes
std::array<Eigen::Vector3d, 3> tri_edge{
{v.col(1) - v.col(0), v.col(2) - v.col(1), v.col(0) - v.col(2)}};
std::array<Eigen::Vector3d, 9> const axx{
{Eigen::Vector3d({0, -tri_edge[0].z(), tri_edge[0].y()}),
Eigen::Vector3d({0, -tri_edge[1].z(), tri_edge[1].y()}),
Eigen::Vector3d({0, -tri_edge[2].z(), tri_edge[2].y()}),
Eigen::Vector3d({tri_edge[0].z(), 0, -tri_edge[0].x()}),
Eigen::Vector3d({tri_edge[1].z(), 0, -tri_edge[1].x()}),
Eigen::Vector3d({tri_edge[2].z(), 0, -tri_edge[2].x()}),
Eigen::Vector3d({-tri_edge[0].y(), tri_edge[0].x(), 0}),
Eigen::Vector3d({-tri_edge[1].y(), tri_edge[1].x(), 0}),
Eigen::Vector3d({-tri_edge[2].y(), tri_edge[2].x(), 0})}};

// Separating axis tests to check if there's a plane separating the
// projections of the AABB and the triangle according to the Separating Axis
// Theorem (see C. Ericson "Real Time Collision Detection" for details)
for (auto const& a : axx)
{
Eigen::Vector3d p = v.transpose() * a;
double const r = e.dot(a.cwiseAbs());
if (std::max(-p.maxCoeff(), p.minCoeff()) > r)
{
return false;
}
}

// Test separating axis corresponding to triangle face normal
Eigen::Vector3d const plane_normal(tri_edge[0].cross(tri_edge[1]));
double const pd = plane_normal.dot(v.row(0));
return testAABBIntersectingPlane(c, e, plane_normal, pd);
}

void markFaults(MeshLib::Mesh& mesh, MeshLib::Mesh const& fault,
int const fault_id, Eigen::Vector3d const& half_cell_size)
{
auto const& elems = mesh.getElements();
std::size_t const n_elems = mesh.getNumberOfElements();
auto mat_ids = mesh.getProperties().getPropertyVector<int>(mat_name);
auto const& fnodes = fault.getNodes();
auto const& felems = fault.getElements();
GeoLib::AABB const fault_aabb(fnodes.cbegin(), fnodes.cend());
auto [min_pnt, max_pnt] = fault_aabb.getMinMaxPoints();

// get bounding box of fault + voxel extent
min_pnt -= half_cell_size;
max_pnt += half_cell_size;

std::array<Eigen::Vector3d, 2> const fault_extent{{min_pnt, max_pnt}};
GeoLib::AABB const fault_aabb_ext(fault_extent.cbegin(),
fault_extent.cend());

// test each voxel grid element vs each fault triangle
Eigen::Vector3d const extent{half_cell_size};
for (std::size_t j = 0; j < n_elems; ++j)
{
// test if bounding box of fault is intersecting voxel
auto const& centre_pnt = MeshLib::getCenterOfGravity(*elems[j]);
if (!fault_aabb_ext.containsPoint(centre_pnt, 0))
{
continue;
}

// test if voxel is intersecting triangle
auto const& c(centre_pnt.asEigenVector3d());
for (auto const* const fault_elem : felems)
{
if (fault_elem->getDimension() != 2)
{
continue;
}

if (testTriangleIntersectingAABB(
*fault_elem->getNode(0), *fault_elem->getNode(1),
*fault_elem->getNode(2), c, extent))
{
(*mat_ids)[j] = fault_id;
break;
}

if (fault_elem->getGeomType() == MeshLib::MeshElemType::QUAD)
{
if (testTriangleIntersectingAABB(
*fault_elem->getNode(0), *fault_elem->getNode(2),
*fault_elem->getNode(3), c, extent))
{
(*mat_ids)[j] = fault_id;
break;
}
}
}
}
}

// test if input mesh is voxel grid
bool isVoxelGrid(MeshLib::Mesh const& mesh)
{
auto const& elements = mesh.getElements();
if (std::any_of(elements.cbegin(), elements.cend(),
[&](auto const& e) {
return (e->getGeomType() !=
MeshLib::MeshElemType::HEXAHEDRON);
}))
{
ERR("Input mesh needs to be voxel grid (i.e. equally sized axis "
"aligned hexahedra).");
return false;
}

for (auto const& e : elements)
{
auto const n = e->getNodes();
if ((*n[0])[2] != (*n[1])[2] || (*n[1])[2] != (*n[2])[2] ||
(*n[4])[2] != (*n[5])[2] || (*n[5])[2] != (*n[6])[2] ||
(*n[1])[0] != (*n[2])[0] || (*n[2])[0] != (*n[5])[0] ||
(*n[0])[0] != (*n[3])[0] || (*n[3])[0] != (*n[7])[0])
{
ERR("Input mesh needs to be voxel grid (i.e. equally sized axis "
"aligned hexahedra).");
return false;
}
}
return true;
}
#include "MeshToolsLib/MeshGenerators/AddFaultToVoxelGrid.h"

int main(int argc, char* argv[])
{
Expand Down Expand Up @@ -236,74 +72,48 @@ int main(int argc, char* argv[])
std::string const fault_name = fault_arg.getValue();
std::string const output_name = output_arg.getValue();

using namespace MeshToolsLib::MeshGenerator::AddFaultToVoxelGrid;

std::unique_ptr<MeshLib::Mesh> mesh(
MeshLib::IO::readMeshFromFile(input_name));
std::unique_ptr<MeshLib::Mesh> fault(
MeshLib::IO::readMeshFromFile(fault_name));
if (mesh == nullptr)
{
ERR("Input mesh not found...");
#ifdef USE_PETSC
MPI_Finalize();
#endif
return EXIT_FAILURE;
}
if (!isVoxelGrid(*mesh))
{
#ifdef USE_PETSC
MPI_Finalize();
#endif
return EXIT_FAILURE;
}
auto const& mat_ids = MeshLib::materialIDs(*mesh);
if (!mat_ids)
{
ERR("Input mesh has no material IDs");
#ifdef USE_PETSC
MPI_Finalize();
#endif
return EXIT_FAILURE;
}

std::unique_ptr<MeshLib::Mesh> fault(
MeshLib::IO::readMeshFromFile(fault_name));
if (mesh == nullptr)
{
ERR("Fault mesh not found...");
#ifdef USE_PETSC
MPI_Finalize();
#endif
return EXIT_FAILURE;
}
if (fault->getDimension() != 2)
{
ERR("Fault needs to be a 2D mesh.");
#ifdef USE_PETSC
MPI_Finalize();
#endif
ERR("Input mesh has no material IDs");
return EXIT_FAILURE;
}

int fault_id = id_arg.getValue();
if (!id_arg.isSet())
{
auto it = std::max_element(mat_ids->cbegin(), mat_ids->cend());
fault_id = *it + 1;
}

Eigen::Vector3d half_cell_size;
if (addFaultToVoxelGrid(mesh.get(), fault.get(), fault_id))
{
auto const n = *mesh->getElement(0)->getNode(0);
auto const c = MeshLib::getCenterOfGravity(*mesh->getElement(0));
half_cell_size[0] = std::abs(c[0] - n[0]);
half_cell_size[1] = std::abs(c[1] - n[1]);
half_cell_size[2] = std::abs(c[2] - n[2]);
MeshLib::IO::VtuInterface vtu(mesh.get());
vtu.writeToFile(output_name);
#ifdef USE_PETSC
MPI_Finalize();
#endif
INFO("The fault was successfully added.");
return EXIT_SUCCESS;
}

markFaults(*mesh, *fault, fault_id, half_cell_size);

MeshLib::IO::VtuInterface vtu(mesh.get());
vtu.writeToFile(output_name);
#ifdef USE_PETSC
MPI_Finalize();
#endif
return EXIT_SUCCESS;
ERR("No fault could be added.");
return EXIT_FAILURE;
}
6 changes: 6 additions & 0 deletions MeshLib/Mesh.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -273,6 +273,12 @@ PropertyVector<int> const* materialIDs(Mesh const& mesh)
return nullptr;
}

PropertyVector<int>* materialIDs(Mesh& mesh)
{
return const_cast<PropertyVector<int>*>(
MeshLib::materialIDs(std::as_const(mesh)));
}

PropertyVector<std::size_t> const* bulkNodeIDs(Mesh const& mesh)
{
auto const& properties = mesh.getProperties();
Expand Down
1 change: 1 addition & 0 deletions MeshLib/Mesh.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -182,6 +182,7 @@ inline bool operator!=(Mesh const& a, Mesh const& b)
/// If the property does not exists (or is of different type), a nullptr is
/// returned.
PropertyVector<int> const* materialIDs(Mesh const& mesh);
PropertyVector<int>* materialIDs(Mesh& mesh);
PropertyVector<std::size_t> const* bulkNodeIDs(Mesh const& mesh);
PropertyVector<std::size_t> const* bulkElementIDs(Mesh const& mesh);

Expand Down
Loading

0 comments on commit 86515bc

Please sign in to comment.