https://gitlab.opengeosys.org/ogs/ogs.git
Tip revision: 10a1660b266d5ce990a32ce7098825b2d898b899 authored by montoyav on 22 March 2021, 16:59:16 UTC
Change of nomenclature (in the documentation) of the physical constants in the Material Library
Change of nomenclature (in the documentation) of the physical constants in the Material Library
Tip revision: 10a1660
SurfaceGrid.cpp
/*
* \date 2012-09-22
* \brief Definition of the SurfaceGrid class.
*
* \file
* \copyright
* Copyright (c) 2012-2021, OpenGeoSys Community (http://www.opengeosys.org)
* Distributed under a Modified BSD License.
* See accompanying file LICENSE.txt or
*/
#include "SurfaceGrid.h"
#include <algorithm>
#include <bitset>
#include "BaseLib/Logging.h"
#include "BaseLib/Error.h"
#include "MathLib/Point3d.h"
#include "Surface.h"
#include "Triangle.h"
namespace GeoLib {
SurfaceGrid::SurfaceGrid(Surface const*const sfc) :
AABB(sfc->getAABB()), _n_steps({{1,1,1}})
{
// enlarge the bounding, such that the points with maximal coordinates
// fits into the grid
for (std::size_t k(0); k<3; ++k) {
_max_pnt[k] += std::abs(_max_pnt[k]) * 1e-6;
if (std::abs(_max_pnt[k]) < std::numeric_limits<double>::epsilon()) {
_max_pnt[k] = (_max_pnt[k] - _min_pnt[k]) * (1.0 + 1e-6);
}
}
std::array<double, 3> delta{{ _max_pnt[0] - _min_pnt[0],
_max_pnt[1] - _min_pnt[1], _max_pnt[2] - _min_pnt[2] }};
if (delta[0] < std::numeric_limits<double>::epsilon()) {
const double max_delta(std::max(delta[1], delta[2]));
_min_pnt[0] -= max_delta * 0.5e-3;
_max_pnt[0] += max_delta * 0.5e-3;
delta[0] = _max_pnt[0] - _min_pnt[0];
}
if (delta[1] < std::numeric_limits<double>::epsilon()) {
const double max_delta(std::max(delta[0], delta[2]));
_min_pnt[1] -= max_delta * 0.5e-3;
_max_pnt[1] += max_delta * 0.5e-3;
delta[1] = _max_pnt[1] - _min_pnt[1];
}
if (delta[2] < std::numeric_limits<double>::epsilon()) {
const double max_delta(std::max(delta[0], delta[1]));
_min_pnt[2] -= max_delta * 0.5e-3;
_max_pnt[2] += max_delta * 0.5e-3;
delta[2] = _max_pnt[2] - _min_pnt[2];
}
const std::size_t n_tris(sfc->getNumberOfTriangles());
const std::size_t n_tris_per_cell(5);
std::bitset<3> dim; // all bits are set to zero.
for (std::size_t k(0); k < 3; ++k)
{
if (std::abs(delta[k]) >= std::numeric_limits<double>::epsilon())
{
dim[k] = true;
}
}
// *** condition: n_tris / n_cells < n_tris_per_cell
// where n_cells = _n_steps[0] * _n_steps[1] * _n_steps[2]
// *** with _n_steps[0] = ceil(pow(n_tris*delta[0]*delta[0]/(n_tris_per_cell*delta[1]*delta[2]), 1/3.)));
// _n_steps[1] = _n_steps[0] * delta[1]/delta[0],
// _n_steps[2] = _n_steps[0] * delta[2]/delta[0]
auto sc_ceil = [](double v){
return static_cast<std::size_t>(std::ceil(v));
};
switch (dim.count()) {
case 3: // 3d case
_n_steps[0] = sc_ceil(std::cbrt(
n_tris*delta[0]*delta[0]/(n_tris_per_cell*delta[1]*delta[2])));
_n_steps[1] = sc_ceil(_n_steps[0] * std::min(delta[1] / delta[0], 100.0));
_n_steps[2] = sc_ceil(_n_steps[0] * std::min(delta[2] / delta[0], 100.0));
break;
case 2: // 2d cases
if (dim[0] && dim[2]) { // 2d case: xz plane, y = const
_n_steps[0] = sc_ceil(std::sqrt(n_tris*delta[0]/(n_tris_per_cell*delta[2])));
_n_steps[2] = sc_ceil(_n_steps[0]*delta[2]/delta[0]);
}
else if (dim[0] && dim[1]) { // 2d case: xy plane, z = const
_n_steps[0] = sc_ceil(std::sqrt(n_tris*delta[0]/(n_tris_per_cell*delta[1])));
_n_steps[1] = sc_ceil(_n_steps[0] * delta[1]/delta[0]);
}
else if (dim[1] && dim[2]) { // 2d case: yz plane, x = const
_n_steps[1] = sc_ceil(std::sqrt(n_tris*delta[1]/(n_tris_per_cell*delta[2])));
_n_steps[2] = sc_ceil(n_tris * delta[2] / (n_tris_per_cell*delta[1]));
}
break;
case 1: // 1d cases
for (std::size_t k(0); k<3; ++k) {
if (dim[k]) {
_n_steps[k] = sc_ceil(static_cast<double>(n_tris)/n_tris_per_cell);
}
}
}
// some frequently used expressions to fill the grid vectors
for (std::size_t k(0); k<3; k++) {
_step_sizes[k] = delta[k] / _n_steps[k];
_inverse_step_sizes[k] = 1.0 / _step_sizes[k];
}
_triangles_in_grid_box.resize(_n_steps[0]*_n_steps[1]*_n_steps[2]);
sortTrianglesInGridCells(sfc);
}
void SurfaceGrid::sortTrianglesInGridCells(Surface const*const sfc)
{
for (std::size_t l(0); l<sfc->getNumberOfTriangles(); l++) {
if (! sortTriangleInGridCells((*sfc)[l])) {
Point const& p0(*((*sfc)[l]->getPoint(0)));
Point const& p1(*((*sfc)[l]->getPoint(1)));
Point const& p2(*((*sfc)[l]->getPoint(2)));
ERR("Sorting triangle {:d} [({:f},{:f},{:f}), ({:f},{:f},{:f}), "
"({:f},{:f},{:f}) into "
"grid.",
l, p0[0], p0[1], p0[2], p1[0], p1[1], p1[2], p2[0], p2[1],
p2[2]);
OGS_FATAL("");
}
}
}
bool SurfaceGrid::sortTriangleInGridCells(Triangle const*const triangle)
{
// compute grid coordinates for each triangle point
std::optional<std::array<std::size_t, 3> const> c_p0(
getGridCellCoordinates(*(triangle->getPoint(0))));
if (!c_p0)
{
return false;
}
std::optional<std::array<std::size_t, 3> const> c_p1(
getGridCellCoordinates(*(triangle->getPoint(1))));
if (!c_p1)
{
return false;
}
std::optional<std::array<std::size_t, 3> const> c_p2(
getGridCellCoordinates(*(triangle->getPoint(2))));
if (!c_p2)
{
return false;
}
// determine interval in grid (grid cells the triangle will be inserted)
std::size_t const i_min(std::min(std::min((*c_p0)[0], (*c_p1)[0]), (*c_p2)[0]));
std::size_t const i_max(std::max(std::max((*c_p0)[0], (*c_p1)[0]), (*c_p2)[0]));
std::size_t const j_min(std::min(std::min((*c_p0)[1], (*c_p1)[1]), (*c_p2)[1]));
std::size_t const j_max(std::max(std::max((*c_p0)[1], (*c_p1)[1]), (*c_p2)[1]));
std::size_t const k_min(std::min(std::min((*c_p0)[2], (*c_p1)[2]), (*c_p2)[2]));
std::size_t const k_max(std::max(std::max((*c_p0)[2], (*c_p1)[2]), (*c_p2)[2]));
const std::size_t n_plane(_n_steps[0]*_n_steps[1]);
// insert the triangle into the grid cells
for (std::size_t i(i_min); i<=i_max; i++) {
for (std::size_t j(j_min); j<=j_max; j++) {
for (std::size_t k(k_min); k<=k_max; k++) {
_triangles_in_grid_box[i+j*_n_steps[0]+k*n_plane].push_back(triangle);
}
}
}
return true;
}
std::optional<std::array<std::size_t, 3>> SurfaceGrid::getGridCellCoordinates(
MathLib::Point3d const& p) const
{
std::array<std::size_t, 3> coords{{
static_cast<std::size_t>((p[0]-_min_pnt[0]) * _inverse_step_sizes[0]),
static_cast<std::size_t>((p[1]-_min_pnt[1]) * _inverse_step_sizes[1]),
static_cast<std::size_t>((p[2]-_min_pnt[2]) * _inverse_step_sizes[2])
}};
if (coords[0]>=_n_steps[0] || coords[1]>=_n_steps[1] || coords[2]>=_n_steps[2]) {
DBUG(
"Computed indices ({:d},{:d},{:d}), max grid cell indices "
"({:d},{:d},{:d})",
coords[0], coords[1], coords[2], _n_steps[0], _n_steps[1],
_n_steps[2]);
return std::optional<std::array<std::size_t, 3>>();
}
return std::optional<std::array<std::size_t, 3>>(coords);
}
bool SurfaceGrid::isPointInSurface(MathLib::Point3d const& pnt,
double eps) const
{
std::optional<std::array<std::size_t, 3>> optional_c(
getGridCellCoordinates(pnt));
if (!optional_c) {
return false;
}
std::array<std::size_t, 3> c(optional_c.value());
std::size_t const grid_cell_idx(c[0]+c[1]*_n_steps[0]+c[2]*_n_steps[0]*_n_steps[1]);
std::vector<Triangle const*> const& triangles(_triangles_in_grid_box[grid_cell_idx]);
auto const it = std::find_if(triangles.begin(), triangles.end(),
[eps, pnt](auto const* triangle) {
return triangle->containsPoint(pnt, eps);
});
return it != triangles.end();
}
} // end namespace GeoLib